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Abstract 

This  thesis  studies  the  problems  of  generation  and  maintenance  of  recirculations  by 
Gulf  Stream  instabilities.  Observations  show  that  the  horizontal  structure  of  the  jet 
and  its  recirculations  suffer  significant  changes  in  time.  Here,  the  role  of  internal 
dynamics  of  the  jet  is  isolated  as  one  of  the  possible  sources  of  such  variability,  and 
the  differences  between  barotropic  and  baroclinic  instabilities  are  investigated. 

The  problem  of  recirculation  development  is  considered  in  a  framework  of  a 
free  spin  down  of  the  2-layer  and  the  l|-layer,  zonally  symmetric,  quasi-geostrophic 
jets.  Linear  stability  analysis  shows  that  in  strongly  baroclinic  basic  flows,  eddies  are 
capable  of  driving  recirculations  in  the  lower  layer  through  the  residual  meridional  cir¬ 
culation.  In  strongly  barotropic  jets,  the  linearly  most  unstable  wave  simply  diffuses 
the  jet.  Nonlinear  stability  analysis  indicates  that  recirculations  are  robust  features 
of  the  2-layer  model.  The  strength  of  recirculations  is  a  function  of  the  model’s  pa¬ 
rameters.  It  increases  with  a  decrease  in  the  value  of  the  nondimensional  0  due  to 
potential  vorticity  homogenization  constrained  by  enstrophy  conservation.  The  recir¬ 
culation  strength  is  a  non-monotonic  function  of  the  baroclinic  velocity  parameter;  it 
is  the  strongest  for  strongly  baroclinic  basic  flows,  weakest  for  flows  with  intermedi¬ 
ate  baroclinic  structure  and  of  medium  strength  for  strongly  barotropic  basic  flows. 
Such  non-monotonic  behavior  is  the  result  of  two  different  processes  responsible  for 
the  recirculation  development:  linear  eddy-mean  flow  interactions  for  strongly  baro¬ 
clinic  basic  flows  and  strongly  nonlinear  eddy-eddy  and  eddy-mean  flow  interaction 
for  strongly  barotropic  flows.  In  the  case  of  the  reduced-gravity  model,  recircula¬ 
tions  develop  only  for  infinite  deformation  raduis.  Basic  flows  with  finite  deformation 
radius  are  only  weakly  supercritical  and  therefore  produced  negligible  recirculations 
after  equilibration. 

The  problem  of  maintenance  of  the  recirculations  is  coupled  to  the  questions  of 
existence  of  low  frequency  variability  and  of  multiple  dynamical  regimes  of  a  system 
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consisting  of  a  quasi-geostrophic  jet  and  its  recirculations.  The  problem  is  studied  in 
a  framework  of  a  2-layer  or  a  reduced-gravity  colliding  jets  model  which  has  no  wind¬ 
forcing.  Instead,  it  is  forced  by  inflows  and  outflows  through  the  open  boundaries. 
Only  the  western  boundary  of  the  domain  is  closed,  and  a  free  slip  boundary  condition 
is  used  there.  The  results  of  the  numerical  experiments  show  that  when  only  the 
mechanism  of  barotropic  instability  is  present,  the  model  has  two  energy  states  for  a 
wide  range  of  interfacial  friction  coefficients.  The  high  energy  state  is  characterized 
by  well-developed  recirculations  and  displays  strong  variability  associated  with  either 
large  recirculating  gyres  and  a  weak  eddy  field  or  small  recirculations  and  a  strong 
eddy  field.  The  low  energy  state  is  characterized  by  large  meridional  excursions  in  the 
separation  point  and  large  amplitude,  westward  propagating  meanders  that  produce 
strong  rings  after  interacting  with  the  western  wall. 

For  physically  relevant  bottom  friction  values,  the  presence  of  baroclinic  in¬ 
stability  in  the  recirculation  regions  of  the  2-layer  model  allows  for  a  unique  dy¬ 
namical  regime  characterized  by  well-developed  recirculations  in  both  layers.  The 
low-frequency  variability  associated  with  the  regime  is  weak  and  is  related  to  merid¬ 
ional  shifts  in  the  position  of  the  jet,  to  wrapping  of  the  recirculations  around  each 
other,  and  to  pulsations  in  their  zonal  extent.  For  strong  bottom  friction,  the  2-layer 
model  has  only  the  mechanism  of  barotropic  instability  which  reduces  it  to  a  l|-layer 
configuration;  the  model  displays  two  dynamical  regimes  and  strong  low  frequency 
variability  in  the  upper  layer,  while  the  lower  layer  is  strongly  frictional. 

Thesis  Supervisor:  Glenn  R.  Flierl 
Title:  Professor 

Massachusetts  Institute  of  Technology 
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Chapter  1 


Introduction 


1.1  Background  and  Motivation 

1.1.1  Observational 

The  discoveries  of  the  intense  recirculating  gyres  south  and  north  of  the  Gulf  Stream 
explained  why  the  transport  of  the  jet  downstream  of  Cape  Hatteras  (150  Sv)  is  five 
times  larger  than  the  maximum  wind-driven  transport  near  the  Straits  of  Florida  (30 
Sv).  The  anticyclonic  recirculation  west  of  50°W  with  the  transport  of  about  60  Sv 
was  first  proposed  by  Worthington  [40].  By  incorporating  both  mass  and  geostrophy 
balance  into  an  inverse  model  applied  to  the  same  data  set,  Wunsch  [41]  and  Wunsch 
and  Grant  [42]  postulated  a  cyclonic  recirculation  west  of  50°W  and  north  of  30°N 
with  the  transport  of  20  Sv. 

By  using  long-term  moored  instruments  near  55°W,  Schmitz  [37]  described  a 
flow  regime  that  consisted  of  narrow  (200  km  wide)  jet-like  currents  with  eastward 
velocities  on  the  flanks  of  the  Gulf  Stream  and  westward  velocities  further  north  and 
south.  These  currents  were  termed  as  “weakly  depth-dependent”  meaning  that  ve¬ 
locity  amplitudes  were  similar  at  thermocline  and  abyssal  depths  and  varied  between 
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Figure  1-1:  Contours  of  zonal  velocity  (cm  s  *) 
son  [32],  ’ 


along  55°  W.  Adapted  from  Richard- 


6  and  10  cm  s  K  Schmitz  noted  that  these  currents  might  be  partly  responsible  for 
enhancing  the  transport  of  the  Gulf  Stream. 

The  combination  of  surface  drifters,  SOFAR  floats  and  current  meters  allowed 
Richardson  [32]  to  further  study  the  vertical  and  horizontal  structure  of  the  Gulf 
Stream  near  55°W.  In  a  region  of  high  eddy  kinetic  energy  and  its  gradient,  he 
found  an  eastward  jet  flanked  by  two  recirculations.  Estimated  at  93  Sv,  the  total 
mean  volume  transport  of  the  Gulf  Stream  was  three  times  larger  than  that  near 
the  Straits  of  Florida.  Such  an  enhancement  was  due  to  additional  41  Sv  and  29 
Sv  carried  respectively  by  the  300  km  wide  northern  and  southern  counter-currents. 
Their  westward  velocities  did  not  change  significantly  with  depth  and  ranged  from  2 
to  10  cm  s-1  (Figure  1-1). 

The  northern  recirculation  was  examined  by  Hogg  et  al.  [14],  Based  on  long¬ 
term  current  and  tracer  measurements,  they  concluded  that  the  Northern  Recirculat¬ 
ing  Gyre,  located  between  the  New  England  Seamount  Chain  and  the  Grand  Banks, 
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Figure  1-2:  Schematic  circulation  in  the  western  North  Atlantic.  Adapted  from  Hogg 
[12]. 

transported  about  40  Sv  of  water  to  the  west  of  63°W.  Half  of  this  water  was  re¬ 
turned  back  to  the  deep  Gulf  Stream.  Consistent  with  the  previous  studies,  the 
vertical  structure  of  the  flow  between  500  m  and  the  bottom  was  weakly  depth  de¬ 
pendent  and  similar  to  that  of  the  southern  recirculation.  It  was  found  also  that  the 
northern  recirculation  exchanged  water  properties  with  the  Deep  Western  Boundary 
Current.  Such  tracers  as  freon  and  oxygen  were  further  advected  into  the  interior  and 
homogenized  within  the  gyre. 

Estimates  of  the  synoptic  transport  of  the  Gulf  Stream  were  obtained  by  Hogg 
[12].  Synoptic  mean  circulation  is  different  from  traditional  Eulerian  mean  circulation 
in  that  it  is  computed  in  a  coordinate  frame  aligned  with  the  instantaneous  axis 
of  the  Stream,  and  therefore  it  accounts  for  changes  both  in  the  path  and  in  the 
meandering  activity  of  the  jet.  Hogg  showed  that  downstream  of  Cape  Hatteras, 
the  synoptic  transport  was  depth-independent  and  reached  150  Sv  between  55°W 
and  60°W.  He  argued  that  about  120  Sv  might  be  evenly  partitioned  between  the 
northern  and  southern  recirculations  in  the  synoptic  mean  frame  compared  to  about 
60  Sv  in  Eulerian  mean  frame.  Figure  1-2  shows  a  schematic  circulation  in  the  western 
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a)  Elongated  Gyre  (December  1986) 


J - 1 - 1 _ , _ L 


70*W  60“W  50"W 


Figure  1-3:  Maps  of  monthly  averaged  sea  surface  height  anomalies  in  western  North 
Atlantic  for  December  1986  (a)  and  October  1988  (b).  Shaded  area  represents  anoma¬ 
lies  larger  than  2.5  m  in  the  southern  recirculating  gyre.  Adapted  from  Kelly  et  al 

[19]- 

North  Atlantic  which  is  characterized  by  the  eastward  flowing  Gulf  Stream  and  two 
recirculations. 

Observations  revealed  that  the  horizontal  structure  of  the  Gulf  Stream  recir¬ 
culations  is  strongly  time-dependent.  For  example,  the  Geosat  altimeter  observations 
of  sea  surface  height  in  the  Gulf  Stream  ( Kelly  et  al.  [19])  have  demonstrated  a 
trend  over  a  period  of  two  years  from  a  state  with  elongated  recirculating  gyres  and 
a  weakly  meandering  jet  to  a  state  with  contracted  gyres  and  a  strongly  meander¬ 
ing  eastward  jet  (Figure  1-3).  There  was  no  significant  correlation  between  the  gyre 
fluctuations  and  the  non-seasonal  wind  curl  line  tilt.  However  the  surface  transport 
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fluctuations,  which  Kelly  et  al.  suggested  were  due  to  variations  in  the  longitudinal 
size  of  the  recirculations,  were  correlated  with  the  changes  in  both  the  jet’s  path  and 
in  the  local  wind  stress  fields.  The  dominant  mode  of  variability  had  a  time  scale  of 
5  to  9  months. 

Observational  data  of  sea  surface  temperature  of  the  Gulf  Stream  between 
75°W  and  60°W  from  April  of  1982  to  December  of  1989  revealed  a  dominant  9- 
month  period  in  the  meandering  intensity  of  the  jet  {Lee  and  Cornillon  [20]).  The 
correlation  between  the  meandering  intensity  of  the  jet  and  the  temporal  variation 
in  the  mean  position  of  the  jet,  which  had  annual  and  inter-annual  oscillations,  was 
small.  Lee  and  Cornillon  suggested  that  the  9-month  cycle  might  be  connected  to 
the  instability  time  scale. 

The  recirculations  might  also  affect  the  dynamics  of  the  Gulf  Stream  by  ad- 
vecting  a  given  meander  westward  and  thus  by  stimulating  the  ring  formation  process. 
On  the  other  hand,  due  to  this  westward  advection,  rings  ultimately  rejoin  the  Gulf 
Stream  and  return  some  of  the  heat  and  potential  vorticity  to  the  jet  {Richardson 
[31]).  Warm  core  rings  which  are  formed  on  the  northern  side  of  the  Stream  bring 
anticyclonic  vorticity  into  the  cyclonic  Northern  Recirculating  gyre.  The  opposite  is 
true  for  the  cold  core  rings  introducing  cyclonic  vorticity  into  the  anticyclonic  south¬ 
ern  recirculating  gyre.  Therefore,  the  existence  of  the  rings  may  generally  tend  to 
decelerate  the  flow  in  the  rectified  regions.  An  increase  in  the  strength  and  amplitude 
of  the  meanders  creates  a  higher  probability  of  having  rings  formed  and,  therefore, 
will  enhance  the  mean  transfer  of  heat  and  momentum  between  the  sub-tropical  and 
the  sub-polar  gyres. 

Low-frequency  variability  is  well  documented  in  observations  of  the  world 
oceans.  Twenty  two  years  of  hydrographic  data  along  137°  show  strong  inter-annual 
bimodal  path  variations  in  the  Kuroshio  {Qiu  and  Joyce  [30]).  The  net  transport  of 
the  Kuroshio  system  which  includes  the  jet  and  its  recirculations  increases  by  30% 
during  the  meander-path  years.  The  increase  is  primarily  due  to  the  decrease  in  the 
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Figure  1-4:  Time  averaged  fields  of  upper  layer  (left)  and  lower  layer  ( middle )  stream- 
function  and  of  interface  height  (right).  Adapted  from  Holland  and  Rhines  [16]. 

recirculation  strength.  Thus,  recirculations  are  strongest  during  the  straight-path 
years.  Observations  of  East  Australian  current  (Roemmich  and  Cornuelle  [33])  and 
of  Brazil  and  Malvinas  currents  (Olson  et  al.  [26])  also  display  strong  inter-annual 
variability. 

1.1.2  Dynamical 

Strong  westward  recirculations  appeared  naturally  in  eddy-resolving,  wind-driven, 
general  circulation  models  (for  a  review  see  Holland  et  al.  [17]).  The  results  of  these 
numerical  experiments  showed  the  importance  of  eddies  in  the  local  potential  vor- 
ticity  balance  of  the  recirculations.  For  example,  in  a  two-layer  quasi-geostrophic 
wind-driven  model  of  Holland  and  Rhines  [16],  two  tight,  intense,  inertial  recircula¬ 
tions  developed  on  the  flanks  of  an  eastward  jet  in  both  layers  (Figure  1-4).  The 
circulation  theorem  applied  to  the  time  mean  streamlines  implied  that  the  abyssal 
recirculations  were  forced  from  above  by  eddy  thickness  (or  eddy  heat)  fluxes,  and 


Figure  1-5:  Time  averaged  upper  layer  (left)  and  lower  layer  ( middle )  eddy  kinetic 
energy  and  eddy  potential  energy  (right).  Adapted  from  Holland  and  Rhines  [16]. 

the  upper  layer  recirculations  were  driven  by  the  eddy  Reynolds  stresses  (or  eddy 
momentum  fluxes).  Energy  transformation  between  the  eddies  and  the  mean  flow 
was  maximum  in  the  areas  of  maximum  eddy  kinetic  and  eddy  potential  energy. 
The  latter  maxima  occurred  in  a  region  of  barotropically  and  baroclinically  unstable 
eastward  jet  and  inertial  recirculations  (Figure  1-5). 

These  results  compared  well  with  observed  surface  eddy  kinetic  (Figure  1- 
6)  and  eddy  potential  energy,  which  were  an  order  of  magnitude  larger  near  strong 
western  boundary  currents  and  their  eastward  extensions  than  in  the  interior  of  the 
ocean.  Holland  and  Rhines  concluded  that  the  abyssal  layer  dynamics  were  primarily 
eddy-driven  and  acted  similar  to  the  down-gradient  diffusion  of  potential  vorticity. 
In  the  upper  layer,  the  same  parameterization  held  true  near  the  jet  separation  point 
and  in  the  region  of  westward  recirculations. 

Some  analytical  models  have  taken  into  account  the  fact  that  the  observed 
transport  of  the  Gulf  Stream  is  several  times  larger  than  that  given  by  the  lin- 
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Figure  1-6:  Surface  eddy  kinetic  energy  per  unit  mass  in  cm2s  2  for  the  North  Atlantic 
Ocean  on  a  1°  grid.  Adapted  from  Schmitz  et  al.  [37]. 

ear  Sverdrup  balance  and  therefore  have  focused  on  strongly  nonlinear  dynamical 
regimes.  The  recirculations  were  given  as  a  steady,  inertial,  almost  free  (forcing  and 
dissipation  were  thought  to  be  negligible  to  the  leading  order)  solution  of  the  quasi- 
geostrophic  equations.  In  both  barotropic  ( Cessi  et  al.  [5])  and  baroclinic  ( Marshall 
and  Nurser  [24];  Greatebatch  [9];  Cessi  [6])  models,  eddies  were  parameterized  as  lat¬ 
eral  down-gradient  diffusion  of  potential  vorticity  and  were  shown  to  play  crucial  role 
in  generating  the  recirculating  gyres. 

Jayne  et  al.  [18]  studied  recirculations  emerging  due  to  both  the  internal  dy¬ 
namics  of  the  jet  and  the  inertial  forcing.  A  quasigeostrophic  barotropic  model  was 
forced  by  an  unstable  eastward  jet  inflow  on  the  western  boundary  and  a  stabilized 
eastward  jet  outflow  on  the  eastern  boundary.  The  resulting  barotropic  eddies  pro¬ 
duced  regions  of  homogenized  potential  vorticity  north  and  south  of  the  jet,  giving 
rise  to  recirculations.  Their  strength  was  found  to  be  a  monotonically  decreasing 
function  of  nondimensional  /?.  However,  it  was  concluded  that  in  the  sense  of  “eddy 
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Sverdrup  relation”,  (3v  =  -V  -uY,  the  recirculations  were  predominantly  inertial 
currents  and  only  weakly  forced  by  the  eddies. 

Therefore,  the  results  of  numerical  experiments  and  the  observational  data 
suggest  strong  correlation  between  the  dynamics  of  the  Gulf  Stream,  its  recirculations 
and  the  western  boundary  current.  In  the  first  part  of  the  thesis,  we  study  the  problem 
of  the  development  of  the  recirculations  during  free  spin  down  of  an  unstable,  zonally 
symmetric,  quasi-geostrophic  jet.  The  following  questions  are  addressed: 

•  Can  recirculations  develop  in  a  zonally  symmetric  unstable  jet?  If  so,  what  is 
the  basic  mechanism  driving  the  recirculations? 

•  How  do  the  jet  strength  and  structure  influence  recirculation  characteristics? 

•  What  are  the  differences  in  the  equilibrated  structure  for  baroclinic  and  barotropic 
jets? 

Zonal  symmetry  combined  with  the  absence  of  any  forcing  excludes  the  possibility 
of  closing  of  the  potential  vorticity  contours  which  may  lead  to  the  recirculation 
development.  Hence,  eddies  generated  by  the  instabilities  if  the  jet  provide  the  only 
mechanism  for  the  changes  in  the  mean  flow  structure  and  therefore  for  the  generation 
of  recirculations. 

Many  time-dependent  wind-driven  models  have  exhibited  intrinsic  low  fre¬ 
quency  variability  that  is  often  associated  with  the  existence  of  multiple  dynamical 
regimes.  For  example,  in  a  reduced-gravity,  free  slip,  quasi-geostrophic  model  with 
interfacial  friction  and  biharmonic  diffusion,  McCalpin  and  Haidvogel  [25]  observed 
chaotic  transitions  between  three  different  regimes  (Figure  1-7):  a  high  energy  state 
with  a  strong,  straight  and  deep  penetrating  eastward  jet,  a  low  energy  state  with  a 
weakly  penetrating  jet  and  violent  meanders  and  eddies,  and  a  medium  energy  state 
with  intermediate  jet  penetration  and  some  eddy/meander  formation.  It  was  hypoth¬ 
esized  that  variations  in  the  eddy-western  boundary  current  interaction  provided  a 
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Medium  Energy  Period 


Figure  1-7:  Contours  of  the  time-averaged  interface  anomaly  field  for  high  energy 
period  (lower  left),  low  energy  period  (lower  right),  and  medium  energy  period  (upper 
right).  Adapted  from  McCalpin  and  Haidvogel  [25]. 

mechanism  for  the  persistence  of  different  regimes  and  the  transitions  between  them. 
Similar  high  and  low  energy  regimes  were  observed  by  Spall  [39]  in  a  wind-forced 
primitive  equation  model  with  bottom  friction  and  lateral  diffusion  and  with  open 
boundary  conditions.  The  presence  of  the  Deep  Western  Boundary  Current  and  its 
interaction  with  the  upper  layer  wind-driven  current  were  crucial  for  the  existence  of 
multiple  dynamical  regimes.  A  number  of  wind-driven,  two-layer,  quasi-geostrophic 
models  also  showed  multiple  dynamical  regimes,  such  as  those  of  Meacham  (personal 
communication)  and  Berio ff  and  McWilliams  [4]. 

However,  in  the  wind-driven  models,  is  it  hard  to  distinguish  between  two 
possible  mechanisms  affecting  low  frequency  variability:  the  global  structure  of  the 
wind  forcing  and  the  internal  dynamics  of  an  unstable  jet.  In  the  present  study, 
we  isolate  the  latter  as  the  mechanism  responsible  for  low  frequency  variability  and 
address  the  following  questions: 
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•  Are  low  frequency  variability  and  multiple  dynamical  regimes  intrinsic  parts  of 
the  internal  dynamics  of  a  quasi-geostrophic  jet  and  its  recirculations? 

•  Are  there  any  differences  between  barotropic  and  baroclinic  instabilities  as 
mechanisms  affecting  low  frequency  variability? 

•  Which  dynamical  regimes  are  associated  with  low  frequency  variability  and 
which  processes  are  responsible  for  the  transitions  between  them? 

•  Are  there  any  differences  between  low  frequency  variability  in  the  globally  forced 
and  boundary-forced  models? 

1.2  Outline  of  the  Thesis 

This  thesis  studying  the  problems  of  generation  and  maintenance  of  recirculations 
by  the  Gulf  Stream  instabilities  consists  of  two  parts.  The  first  part  investigates  the 
development  of  recirculations  during  the  spin  down  of  a  quasi-geostrophic,  zonally 
symmetric,  unstable  jet.  The  second  part  of  the  thesis  examines  the  problem  of 
existence  and  maintenance  of  low  frequency  variability  and  of  multiple  dynamical 
regimes  due  to  the  mechanisms  of  barotropic  and  baroclinic  instabilities. 

Chapter  2  addresses  the  first  problem  in  a  framework  of  the  reduced-gravity 
and  the  two-layer  quasi-geostrophic,  doubly-periodic,  zonally  symmetric  models.  It 
examines  the  following  questions:  what  is  the  basic  mechanism  for  driving  the  west¬ 
ward  flows;  how  do  the  jet  strength  and  structure  influence  recirculation  character¬ 
istics;  what  are  the  differences  in  the  equilibrated  structure  for  the  baroclinic  and 
barotropic  jets?  The  chapter  begins  with  the  formulation  of  the  numerical  model 
and  is  followed  by  the  discussion  of  the  results  of  a  time-dependent  two-layer  case. 
The  strength  and  the  structure  of  the  recirculations,  arising  due  to  eddy  heat  and 
momentum  fluxes,  are  examined  for  different  values  of  nondimensional  (3— parameter 
and  for  a  set  of  basic  flows  with  varying  degree  of  baroclinicity.  To  study  the  nature  of 
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eddy-mean  flow  interaction  and  its  role  in  recirculation  generation,  a  linear  stability 
analysis  is  performed  for  basic  flows  with  different  degrees  of  baroclinicity  and  various 
jet  widths.  A  hierarchy  of  “mean-flow”  models,  which  includes  only  a  mechanism  of 
eddy-mean  flow  interactions,  is  further  examined.  The  formulation  and  the  numeri¬ 
cal  results  of  the  time-dependent  reduced-gravity  model  and  a  corresponding  linear 
stability  analysis  follow  the  discussion. 

Chapters  3  introduces  a  two-layer,  quasi-geostrophic,  colliding  jets  model, 
which  is  forced  by  the  inflows  and  outflows  through  the  open  boundaries.  Such  forc¬ 
ing  allows  the  isolation  of  barotropic  and  baroclinic  instabilities  as  the  mechanisms 
responsible  for  multiple  dynamical  states  and  low  frequency  variability.  The  chapter 
begins  with  the  description  of  model  equations,  initial  and  boundary  conditions,  fric¬ 
tional  operators  and  dimensional  and  nondimensional  parameters.  The  importance  of 
local  forcing  and  its  difference  from  the  global  wind  forcing  are  illustrated  in  a  section 
on  model  energetics.  The  dynamics  of  rings  on  the  western  boundary  are  discussed 
at  the  end  of  the  chapter. 

Chapter  4  reports  the  results  of  the  reduced-gravity,  colliding  jets  model,  which 
isolates  the  importance  of  barotropic  instability  and  dissipation  for  low  frequency  vari¬ 
ability.  Following  the  model  formulation,  the  antisymmetric  steady  solutions  and  their 
linear  stability  properties  are  presented  for  various  values  of  the  interfacial  drag  coef¬ 
ficient.  The  chapter  further  focuses  on  the  sensitivity  of  the  time-dependent  solutions 
to  changes  in  biharmonic  diffusion  and  interfacial  friction.  A  detailed  description  of 
a  reference  run  shows  that  the  model  has  two  preferred  energy  states,  each  associ¬ 
ated  with  a  unique  dynamical  regime.  The  discussion  of  the  dynamical  differences 
between  the  states  and  their  transitions  as  well  as  a  comparison  with  the  results  of 
the  wind-forced  models  follow. 

Chapter  5  addresses  a  problem  of  low  frequency  variability  and  multiple  dy¬ 
namical  states  in  a  framework  of  a  two-layer,  colliding  jets  model,  which  includes 
both  baroclinic  and  barotropic  instability  mechanisms.  The  sensitivity  of  the  model 
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to  changes  in  the  value  of  the  bottom  friction  parameter  is  examined  and  is  related 
to  a  scaling  analysis.  The  discussion  of  a  reference  run  includes  the  analysis  of  energy 
time  series,  an  extended  empirical  orthogonal  function  decomposition  for  low-  and 
high-frequency  variability  and  a  diagnosis  of  zones  of  baroclinic  growth.  The  results 
are  compared  with  those  of  the  reduced-gravity,  colliding  jets  model  and  a  two-layer, 
wind-forced  models. 

Chapter  6  presents  conclusions  for  the  Thesis,  summarizes  limitations  of  the 
present  study,  poses  questions  for  further  research  and  discusses  relevance  of  this  work 
to  the  real  ocean. 


23 


Chapter  2 


Generation  of  Recirculations  in  a 
Two-Layer  Quasi-Geostrophic 
Zonally  Symmetric  Jet  Model 


The  present  Chapter  addresses  the  following  questions: 

•  Can  recirculations  develop  in  a  zonally  symmetric  unstable  jet?  If  so,  what  is 
the  basic  mechanism  driving  the  recirculations? 

•  How  do  the  jet  strength  and  structure  influence  recirculation  characteristics? 

•  What  are  the  differences  in  the  equilibrated  structure  for  baroclinic  and  barotropic 
jets? 
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2.1  Model  Formulation 


2.1.1  Model  Equations 


A  time-dependent,  two-layer,  quasi-geostrophic,  semi-spectral  model  is  integrated  in 
a  square,  doubly  periodic  domain.  The  governing  potential  vorticity  equations  are 


,d_  dpi  d 
dt  dx  dy 


dpi  d 
dy  dx 


)qi  +  P 


dpi 

dx 


—  F  +  <Sj, 


(2.1) 


where  q{  =  V2pt  -  Fi(pi  -  p^-i),  F{  =  f^L2/g'Du  j3  =  pdimL2  jU  and  i  =  1, 2.  Since 
we  are  interested  in  the  problem  of  a  nonlinear  equilibration  of  a  quasi-geostrophic 
jet,  the  forcing  is  set  to  zero.  However,  a  numerical  filter,  T,  is  required  for  numerical 
stability;  it  acts  to  remove  small  scale  enstrophy  without  dissipating  too  much  energy. 
Because  of  the  periodicity  in  the  meridional  direction,  two  sponge  layers  are  set  up 
near  the  northern  and  southern  boundaries.  In  Equation  (2.1),  they  are  represented 
by  a  sponge  operator  S{.  Sponges  prevent  any  disturbance  which  leaves  the  domain 
through  the  northern  (southern)  boundary  from  reentering  it  from  the  south  (north). 
More  details  on  the  form  of  the  sponge  operator  are  presented  later.  The  numerical 
code  is  modified  from  a  program  by  Glenn  Flierl  (personal  communication). 


2.1.2  Dimensional  Scales 

The  dimensional  scales  of  the  model  are  chosen  as  follows.  The  horizontal  length 
scale,  £,  is  407  km.  The  model  domain  is  square  defined  by  Ldomain  x  Ldomain ,  where 
Ldomain  =  27tL  =  2560  km.  The  upper  and  lower  layer  depths  are  Dx  —  1  km  and 
£>2=4  km  respectively.  The  baroclinic  deformation  radius,  Ldef,  is  40  km.  The 
velocity  scale,  U  —  0.8  m  sec-1,  is  chosen  as  the  maximum  velocity  in  the  upper  layer. 
The  time  scale  is  chosen  as  an  advective  time  scale,  Tadv  —  L/JJ  —  29.6  days.  The 
Coriolis  parameter  is  represented  by  /  =  /o  +  f3y,  where  /0  =  10-4  sec-1  and 
/3  —  2  •  10-u  m-1sec-1. 
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2.1.3  Nondimensional  Parameters 


The  model  is  characterized  by  four  nondimensional  parameters: 

•  the  ratio  of  the  advective  time  scale  to  the  barotropic  Rossby  wave  period: 

P  =  PdimL2/U- 

•  the  inverse  deformation  radius,  defined  as  the  ratio  of  a  horizontal 

length  scale  to  the  internal  baroclinic  deformation  radius:  7  =  \fF[  +  F2  = 

I'/l'de/! 

•  the  depth  ratio:  S  =  Di/D2, 

•  the  ratio  of  the  upper  layer  initial  velocity  amplitude,  uu  to  the  lower  layer 
initial  velocity  amplitude,  u2  :  a  =  Ui/u2. 

In  all  the  experiments,  the  values  of  8  and  7  are  fixed  at  0.25  and  10.18 
respectively  and  the  values  of  /?  and  a  are  varied. 

2.1.4  Sponge  Operator 

As  was  mentioned  before,  two  sponge  layers  are  set  up  near  the  northern  and  southern 
boundaries  to  remove  periodicity  in  the  meridional  direction.  After  performing  a  series 
of  experiments  with  different  widths  of  the  sponge  layer,  we  found  the  optimal  width 
of  about  300  km,  which  allows  enough  distance  for  the  decay  of  the  disturbances. 
Whenever  a  sponge  layer  was  too  narrow,  it  acted  as  a  reflector  rather  than  an 
absorber. 

The  sponge  operator  has  the  following  form: 

Si  =  -nv  (2.2) 

The  first  term  is  represented  by  a  sponge  function,  H  (Figure  2-1),  acting  to  relax 
the  potential  vorticity  field  to  the  planetary  vorticity.  This  type  of  damping  is  similar 
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Figure  2-1:  Sponge  function,  TZ,  defines  the  ratio  of  the  advective  time  scale  to  the 
decay  time  scale.  The  latter  varies  from  being  infinite  in  the  interior  of  the  domain 
to  0.5  day  near  the  boundary. 


to  bottom  friction.  The  meridionally  dependent  friction  coefficient,  7 Z(y),  is  the  ratio 
of  the  advective  time  scale  to  the  decay  time  scale,  Tadv/Tdecay.  Within  each  sponge 
layer,  the  decay  time  scale,  Tdecay ,  slowly  decreases  from  being  infinite  to  0.5  day  at 
the  boundary. 

The  second  term  in  the  sponge  operator  involves  the  meridional  gradient  of 
a  sponge  function  acting  on  the  velocity  field.  It  helps  to  remove  spurious  energy 
arising  from  boundary  currents  with  zero  vorticity.  This  term  comes  naturally  from 
the  following  velocity  damping  equations: 

^  =  ~n(x,  y)(ui  -  U 0,  *  =  1,  2,  (2.3) 


where  u*  =  {ui,Vi)T  is  a  horizontal  velocity  vector,  U*  =  ( Ui,Vj)T  is  a  relaxation 
velocity  vector  and  TZ  is  a  friction  coefficient,  which  in  general  could  be  a  function  of 
spatial  variables.  The  equation  for  relative  vorticity,  G  —  ^  ~  becomes 
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Zonal  symmetry  in  the  model  implies  that  the  frictional  coefficient  changes  only  in 
the  meridional  direction,  i.e.  =  0.  The  relaxation  velocity,  Ui,  was  chosen  to  be 
zero. 


2.1.5  Basic  Flow 


The  basic  velocity  profile  corresponds  to  a  thin  eastward  jet  and  is  defined  as 

d^(y)  m 


Ui(y )  = 


-,  i  =  1,2. 


(2.5) 


dy  cosh  2wy ' 

The  amplitude  of  the  lower  layer  nondimensional  velocity,  w2,  varies  from  being  zero 
to  the  amplitude  of  the  upper  layer  nondimensional  velocity,  ux.  Thus,  we  can  create 
both  strongly  baroclinic  jets  or  purely  barotropic  jets  (Figure  2-2a).  The  parameter 
w  defines  the  width  of  the  jet  and  is  chosen  to  equal  8.14,  which  corresponds  to  a  100 
km  wide  basic  flow. 

The  basic  potential  vorticity  is  then  given  by 

Qi  =  =  l3v+tmhwy{^^/  +  ~ (“j-*) }  <2-6) 

The  basic  potential  vorticity  gradient  changes  sign,  so  that  it  satisfies  the  necessary 
condition  for  instability  (Figure  2-2b).  The  results  of  the  linear  stability  analysis  are 
presented  later. 


2.1.6  Initialization 

We  perturb  the  basic  flow  given  by  Equation  (2.5)  with  a  sum  of  wave  packets  of 
different  wave  lengths,  each  having  a  Gaussian  envelope.  Initially,  the  disturbance  is 
localized  in  the  middle  of  the  domain. 


2.2  Linear  Stability  Analysis  of  a  Two-Layer  Model 

Two  processes  play  major  roles  in  a  nonlinear  spin  down  of  an  unstable  quasi- 
geostrophic  jet:  eddy-mean  flow  and  eddy-eddy  interactions.  The  following  section 
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Figure  2-2:  The  nondimensional  velocity  (a)  and  potential  vorticity  (b)  for  the  upper 
layer  (solid  line)  and  lower  layer  (dash  line)  initial  basic  flow.  Parameters:  U\  — 
1,  u2  =  0.1,  w  =  8.14,  (5  =  4.15,  7  =  10.18,  S  =  0.25,  U  =  0.8  m  sec-1.  North-south 
extent  of  the  full  domain  is  [— 7T,  7 r]. 

reports  the  results  of  a  linear  stability  analysis  which  studies  the  influence  of  the 
barotropic  and  baroclinic  eddy  field  on  the  mean  flow  structure. 

Although  the  nonlinear  model  defined  by  Equation  (2.1)  has  a  numerical  fric¬ 
tional  operator  and  sponge  layers  near  the  doubly-periodic  meridional  boundaries,  the 
linearized  model  is  solved  for  an  inviscid  case  in  a  channel.  Such  simplifications  are 
possible  for  two  reasons.  First,  in  the  nonlinear  model,  friction  is  very  weak.  Second, 
as  we  show  later,  unstable  normal  modes  resulting  from  the  linear  stability  analysis 
are  zero  in  the  regions  where  a  sponge  operator  of  the  nonlinear  model  is  not  trivial. 
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2.2.1  Linear  Model  Formulation 

Linearizing  the  inviscid  quasi-geostrophic  equations  of  motion  (2.1)  around  a  basic 
flow  gives: 

r\ 

jj-  +  +  J(fa,Qi)  =  0,  (2.7) 

where  %(y)  =  -  Jy  Ui(y')dy'  and  Qi(y)  =  (3y  +  3_,  -  %)  +  are  the  basic  flow 

streamfunction  and  potential  vorticity  fields,  fa  =  fa(x,y,t)  and  q{  =  qi(x,y,t)  are 
the  perturbation  streamfunction  and  potential  vorticity  fields  and  i  =  1,2  is  a  layer 
index.  We  look  for  a  zonally  periodic  solution  of  the  form 

fa(x,y,t)=  Real ($*(?/)  eifc(l_ci)),  (2.8) 

where  A;  is  a  horizontal  wave  number  and  c  is  the  phase  speed  of  the  disturbance. 
Substituting  such  a  solution  into  Equation  (2.7)  gives  the  following  problem  for  the 
normal  mode  meridional  structure: 

(Ui(y)  -  C)(±±  -  k2<$>i  +  Fi($ a-,-  -  $i))  +  =  0,  *  =  1,  2,  (2.9) 

subject  to  no-normal-flow  boundary  conditions  on  the  lateral  boundaries 

$*  (y  =  tt)  =  (y  =  -7 r)  =  0.  (2.10) 

2.2.2  Method  of  Solution 

The  linear  problem  (2.9)  is  solved  numerically.  Discretizing  the  equations  at  yj  = 
-7r  +  Ay(j  -  1),  where  A y  =  ~~  and  j  =  1  , iV,  using  a  second  order  finite  dif¬ 
ference  scheme  and  incorporating  boundary  conditions  (2.10),  provides  a  generalized 
eigenvalue  problem  of  the  form: 

,4$  =  CB$.  (2.11) 

Matrices  A  and  B  have  the  following  structures: 


B=[*  V 

Z,  P2  _ 

where  %i  and  Vt  are  tri-diagonal  matrices  defined  as 

hf]  U[2)  0  0  •••  0  0 

u\z)  hf]  uf]  0  •••0  0 

0  U, }4)  h(4)  u\4)  •••  0  0 

nr  = 

o  o  o  o  •••  hrz>  [/f-2) 

0  0  0  0  •••  E^_1) 

with  diagonal  elements 

hf  =  >  (-2  -  Ay2(fc2  +  F<))  +  Ay2^“ 

U F>  =  f/j(y  =  y2), 

and 

pj2)  1  0  0  •••  0 

1  p{3)  1  0  •••  0 

0  1  pj4)  1  •••0 

0  0  0  0  pjN~2) 

0  0  0  0  •••  1 

with  diagonal  elements 

PiJ)  =  -2  -  A y2(k2  +  Fi) 

for  i  =  1,2  and  j  =  2,  •  •  • ,  N  —  1.  Matrices  Qi  and  are  diagonal  with  nonzero 
elements  equal  g-J)  =  U-^FiAy2  and  =  F{Ay 2  respectively. 

MATLAB’s  SPTARN  routine  was  used  to  solve  the  generalized  eigenvalue 
problem  defined  in  Equation  (2.11).  We  performed  a  series  of  experiments  with 
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RESOLUTION,  N  RESOLUTION,  N 

Figure  2-3:  Real,  cr,  (a)  and  imaginary,  a,  (b)  parts  of  an  eigenvalue  as  functions  of 
resolution,  N,  for  ut  =  1,  u2  =  0.99,  /?  =  4.15,  7  =  10.18,  5  =  0.25,  k  =  8,  U  =  0.8 
m  sec-1. 

different  resolutions  to  ensure  sufficient  convergence.  Figure  2-3  shows  the  real  and 
imaginary  parts  of  an  eigenvalue  as  functions  of  resolution,  N,  when  ux  =  1,  u2  = 
0.99,  (3  =  4.15,  7  =  10.18,  and  k  =  8.  N  —  500  was  found  to  be  sufficient  for  the 
choice  of  parameters  considered  in  the  numerical  experiments. 

2.2.3  Maximum  Growth  Rates  of  Unstable  Modes 

The  first  set  of  experiments  was  performed  for  the  basic  flow  with  the  following 
parameters:  U\  =  1,  (3  =  4.15,  7  =  10.18,  5  =  0.25  and  the  lower  layer  basic  velocity 
amplitude  u2  varying  between  0  and  0.99.  Typical  unstable  eigenvalues  are  presented 
in  Figure  2-4.  As  k  increases,  the  number  of  unstable  modes  decreases  from  three  to 
zero.  Long  waves  are  weakly  unstable.  Notice  that  the  linear  stability  analysis  was 
performed  only  for  the  integer  horizontal  wave  numbers,  k  =  1,2,  ...,n.  The  reason 
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is  that  only  waves  with  such  wave  numbers  are  present  in  the  corresponding  fully 
nonlinear  doubly-periodic  model. 


0  2  4  6  8  10  12  14  16  18  20 

k 


Figure  2-4:  The  growth  rates,  kci,  (o)  and  the  corresponding  real  parts  of  the  phase 
speed,  cr,  (*)  as  functions  of  a  horizontal  wave  number,  k.  Parameters:  ui  =  1,  u2  = 
0.99,  /3  =  4.15,  7  =  10.18,  8  =  0.25,  U  =  0.8  m  sec”1. 


Figure  2-5  shows  the  maximum  growth  rate,  k^,  as  a  function  of  a  horizontal 
wave  number,  k.  In  this  experiment,  the  amplitude  of  the  upper  layer  nondimen- 
sional  velocity,  u\,  is  fixed  at  1.0.  The  amplitude  of  the  lower  layer  nondimensional 
velocity,  u2,  varies  between  0.0  and  0.99.  Clearly,  the  maximum  growth  rate  is  not  a 
monotonic  function  of  u2  (Figure  2-6);  as  the  flow  becomes  more  barotropic,  it  first 
decreases,  reaching  a  minimum  at  u2  =  0.2,  and  then  increases,  reaching  a  maxi¬ 
mum  at  u2  =  0.99.  This  non-monotonic  behavior  suggests  that  it  is  the  baroclinic 
mode  that  dominates  linear  growth  for  u2  <  0.2.  As  the  barotropicity  of  the  basic 
flow  increases,  the  barotropic  mode  becomes  more  and  more  important,  and  finally 
it  dominates  the  linear  growth. 

The  dimensional  length  of  the  most  unstable  wave  increases  from  320  km 
(k  =  6)  to  426  km  ( k  =  8)  with  a  decrease  in  baroclinicity.  Correspondingly,  the 
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(a)  (b) 


Figure  2-5:  The  maximum  growth  rate,  kcu  as  a  function  of  a  horizontal  wave  number, 
k.  The  amplitude  of  the  upper  layer  nondimensional  basic  velocity,  u\,  is  1.  The 
amplitude  of  the  lower  layer  nondimensional  basic  velocity,  u2,  changes  from  0  to  0.2 
(a)  and  from  0.2  to  0.99  (b).  Parameters:  (3  =  4.15,  7  =  10.18,  6  =  0.25,  U  =  0.8m 
sec-1. 

dimensional  phase  speed  changes  from  0.07  m  sec-1  to  0.34  m  sec-1.  Independent  of 
baroclinicity,  the  flow  is  unstable  for  0  <  k  <  15,  so  that  only  waves  shorter  than  170 
km  are  stable. 

The  meridional  structure  of  the  most  unstable  mode  is  similar  for  all  considered 
values  of  the  lower  layer  nondimensional  velocity  amplitude,  u2  (Figure  2-7).  The 
mode  is  nonzero  only  in  the  interior  of  the  domain,  which  justifies  the  use  of  the 
linear  model  in  a  periodic  channel  and  neglecting  of  the  sponge  layers. 

Nondimensional  (3  has  a  stabilizing  effect  on  the  mean  flow  (Figure  2-8).  For 
strongly  baroclinic  cases  (a),  waves  longer  that  900  km  are  linearly  stable  when  (3  = 
6.64.  The  long  wave  cutoff  disappears  for  smaller  (3.  For  barotropic  basic  flows  (b), 
there  is  no  long  wave  cutoff.  Waves  shorter  than  160  km  are  stable. 
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Figure  2-6:  The  maximum  growth  rate,  k,C{,  (solid  line)  and  the  corresponding  nondi- 
mensional  phase  speed,  cr,  (dashed  line)  as  functions  of  the  amplitude  of  the  lower 
layer  nondimensional  basic  velocity,  «2-  Parameters:  U\  =  1.0,  (3  =  4.15,  7  = 
10.18,  8  =  0.25,  U  =  0.8  m  sec-1. 

2.2.4  Wave  -  Zonal  Mean  Flow  Interactions 

The  Eulerian  zonal  mean  circulation  is  calculated  to  understand  how  the  most  linearly 
unstable  mode  affects  the  mean  flow.  We  assume  the  amplitude  of  the  linear  mode 
to  be  small  and  take  the  zonal  average  of  Equations  (2.7).  Repeating  the  arguments 
given  by  Shepherd  [35]  and  adopting  his  notations  imply  the  following  problem: 

+  F(?s-i  -  A))  =  i  =  1,2,  (2.14) 

subject  to  boundary  conditions: 

ipi  =  0  at  y  =  — 7r,7 r,  (2.15) 

where  bar  denotes  zonally  averaged  fields,  and  0  and  q  are  the  streamfunction  and 
the  potential  vorticity  fields  of  the  most  unstable  mode. 
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Figure  2-7:  The  real  (solid  line)  and  imaginary  (dashed  line)  parts  of  the  most  unstable 
mode,  $(y,  k  =  6)  in  the  upper  layer  (a)  and  in  the  lower  layer  (b).  Parameters: 
ui  =  1.0,  u2  =  0.0,  0  =  4.15,  7  =  10.18,  5  =  0.25,  U  =  0.8  m  sec"1. 

Substituting  a  normal  mode  form  into  the  forcing  term,  J(<f>hqi),  of  the  above 
equations  implies  that  the  solution  must  obey  the  following: 

— f/V  =  2kcie(y)e2t«l.  (2.16) 

Thus,  the  problem  (2.14)  can  be  rewritten  as 

ewv  +  Wa-i  -  ©i)  =  —^Im{<j>*iy<t>iyy  -  (j)lyyy(j>t  +  Fi(<f> 3_i(f>*y  -  ^3„%^)),(2.17) 

where  *  =  1,2,  j  =  3  -  *  are  layer  indices,  k  and  kc{  are  the  horizontal  wave  number 
and  the  growth  rate  of  the  most  unstable  mode  and  a  star  sign  denotes  complex 
conjugate.  The  boundary  conditions  are  given  by  0(-7r)  =  0(7r)  =  0. 

Finite  differencing  reduces  problem  (2.17)  to  solving  AO  =  B,  where 

„  & 

A  = 

Qi  7^2 
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(a) 


Figure  2-8:  The  maximum  growth  rate,  kci,  as  a  function  of  a  horizontal  wave  number, 
k:  for  p  =  (6.64,4.74,3.68,2.76,2.21).  Parameters:  U\  =  1,  u2  =  0  in  (a)  and  u\  — 
1,  u2  =  0.99  in  (b). 
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Individual  blocks  of  the  matrix  have  the  following  structures: 


where  ht-  =  -  2  -  A y2Fi,  and 


Qi  = 


where  g{  =  A y2Fi:  i  =  1,2. 


4  -1  0  •••  0  0  0 

hi  1  0  •••  0  0  0 

1  hi  1  •••  0  0  0 
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2.2.5  Eddy  Energy  Equation 


Transient  eddies  can  grow  at  the  expense  of  kinetic  or  potential  energy,  or  both, 
which  implies  barotropic,  baroclinic  or  mixed  nature  of  the  instability  process.  The 
equation  for  eddy  energy,  E  ,  ( Pedlosky  [27])  for  the  linear  model  is  given  by 

~E'  =  J\eKECx  +  EKEC2  +  EPEC )  dy,  (2.18) 


where  the  i— the  layer  eddy  kinetic  energy  conversion,  EKECi,  is  produced  by  the 
depth-weighted  Reynolds  stresses  acting  on  the  basic  state  shear: 


EKECi  = 


Pi  dfa  d<f>j  dUj 
D  dx  dy  dy  ’ 


*  =  1,2, 


(2.19) 
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Figure  2-9:  The  eddy  energy  convergence  terms:  upper  layer  kinetic,  EKEC\,  lower 
layer  kinetic,  EKEC2,  potential,  EPEC ,  and  total  for  u2  =  0  (a),  u2  =  0.2  (b)  and 
if2  =  0.99  (c).  Parameters:  u\  =  1,  (5  —  4.15,  7  =  10.18,  5  =  0.25,  U  —  0.8  m  sec-1. 


and  the  eddy  potential  energy  conversion,  EPEC ,  is  produced  by  the  meridional  eddy 
heat  flux  multiplied  by  the  basic  state  temperature  gradient: 


EPEC 


F\D\  +  F2D2 

2D 


(Ui 


TT  \d(t>2A, 


(2.20) 


When  m2  =  0.99  (Figure  2-9c),  the  instability  is  purely  barotropic;  the  potential 
energy  conversion  is  negligible.  For  u2  =  0.2  (Figure  2-9b),  the  instability  is  mixed; 
both  the  kinetic  and  potential  eddy  energy  conversion  terms  are  positive  and  both 
participate  equally  in  the  growth  in  time  of  the  total  eddy  energy.  For  the  case  of 
«2  =  0  (Figure  2-9a),  the  instability  is  predominantly  baroclinic  with  some  barotropic 
contribution  in  the  upper  layer.  In  the  lower  layer,  the  eddy  energy  increases  only 
due  to  the  eddy  heat  fluxes.  Thus,  taking  into  account  different  nature  of  linear 
instability,  we  expect  different  impact  of  transient  disturbances  on  the  mean  flow. 


2.2.6  Transformed  Eulerian  Mean  Circulation 

To  illuminate  the  dynamics  of  the  linearly  growing  disturbance  and  its  impact  on  the 
mean  flow,  the  Transformed  Eulerian  Mean  (TEM)  circulation  is  calculated  ( Andrews 
and  McIntyre  [1],  Pedlosky  [27],  Shepherd  [35]): 

^  =  v*  +  V'EPj,  i  =  l,2.  (2.21) 

In  the  TEM  formulation,  the  time  variations  of  the  mean  flow  velocity  result 
from  the  divergence  of  the  Eliassen-Palm  fluxes: 

Q 

V-EPj  =  Wi  =  -t^(uFh)  +  (-1  YFiVi^i  -  (j)2 ),  (2.22) 

and  from  the  residual  meridional  circulation: 

vf  =  vf~  (-lyFiVtifa  -cl)2),  (2.23) 

where  vf  is  an  ageostrophic  meridional  velocity  and  ( Ui,Vi)T  =  (~ ^f)r  is  a 
geostrophic  perturbation  velocity  vector. 
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Figure  2-11:  (a):  Upper  layer  zonal  mean  acceleration,  (c):  Upper  layer  north¬ 
ward  potential  vorticity  flux,  W[q{.  (e):  Upper  layer  residual  meridional  circulation, 
vi-  (b),  (d)  and  (f):  same  for  the  lower  layer.  Parameters:  u\  —  1,  u2  =  0.2,  /3  = 
4.15,  7  =  10.18,  S  =  0.25,  U  =  0.8  m  sec-1. 


(a)  (b) 


Figure  2-12:  (a):  Upper  layer  zonal  mean  acceleration,  (c):  Upper  layer  north¬ 
ward  potential  vorticity  flux,  W[q f.  Parameters:  U\  =  1,  u2  =  0.99,  /?  =  4.15,  7  = 
10.18,  S  =  0.25,  U  —  0.8  m  sec-1. 

Figure  2-10a  shows  mean  zonal  accelerations,  ^  =  —2 kCiQiye2kCit ,  for  the 
case  when  the  amplitudes  of  the  upper  and  lower  layer  basic  velocities  are  1.0  and 
0.0  respectively.  Clearly,  the  most  unstable  linear  wave  has  a  “diffusive”  effect  on 
the  upper  layer  mean  flow;  it  broadens  the  flow  by  decreasing  its  zonal  velocity  in 
the  middle  of  the  domain  and  increasing  it  on  the  flanks  of  the  jet.  In  the  lower 
layer,  the  linearly  growing  wave  tends  to  induce  an  eastward  jet  in  the  middle  of  the 
domain  surrounded  by  weak  westward  flows  on  the  north  and  the  south.  Notice  that 
the  amplitude  of  the  lower  layer  zonal  mean  acceleration  is  an  order  of  magnitude 
smaller  than  the  amplitude  in  the  upper  layer. 

This  contrasts  with  the  mean  zonal  accelerations  for  stronger  barotropic  flows. 
For  both  u2  =  0.2  (Figure  2-lla)  and  u2  =  0.99  (Figure  2-12  a),  eddies  have  purely 
diffusive  effect  on  the  mean  flow  in  both  layers.  There  is  no  tendency  for  recirculations 
in  these  cases. 
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As  Equation  (2.21)  shows,  the  mean  zonal  acceleration  is  defined  by  the  inter¬ 
play  between  the  residual  meridional  circulation  and  the  northward  eddy  potential 
vorticity  fluxes.  When  u2  =  0,  they  have  opposite  impact  on  the  mean  flow  in  the  cen¬ 
ter  of  the  jet;  v'q  tends  to  decelerate  the  mean  flow  in  the  upper  layer  and  accelerate 
it  in  the  lower  layer,  whereas  v*  tends  to  decrease  the  isopycnal  slope,  and  therefore 
it  is  northward  (accelerating  the  mean  flow)  in  the  upper  layer  and  southward  (decel¬ 
erating  the  mean  flow)  in  the  lower  layer.  However,  the  dominant  contribution  comes 
from  the  eddy  potential  vorticity  fluxes.  On  the  flanks  of  the  jet  though,  both  Vcf 
and  v*  tend  to  induce  an  eastward  flow  in  the  upper  layer,  so  that  they  both  have 
diffusive  effect  of  the  mean  flow.  In  the  lower  layer,  their  role  is  different:  it  is  only 
due  to  the  residual  meridional  circulation  that  there  is  a  tendency  for  westward  flows 
on  the  flanks  of  the  jet  in  the  lower  layer.  As  the  amplitude  of  the  lower  layer  nondi- 
mensional  basic  velocity  increases,  the  slope  of  the  isopycnal  diminishes.  The  latter 
leads  to  a  substantial  decrease  in  the  strength  of  the  residual  meridional  circulation, 
so  that  it  cannot  dominate  the  eddy  potential  vorticity  fluxes.  Therefore,  there  is  no 
tendency  for  the  formation  of  recirculations  in  the  lower  layer. 

For  strongly  barotropic  basic  flows,  the  residual  meridional  circulation  becomes 
negligible  and  the  eddy  potential  vorticity  fluxes  have  a  purely  diffusive  effect  on  the 
mean  flow  during  the  stage  of  linear  growth. 

2.2.7  Relevance  to  Fully  Nonlinear  Calculations 

As  will  be  described  in  section  2.3,  the  upper  layer  mean  zonal  flow  decreases  its  am¬ 
plitude  and  becomes  broader  during  the  first  two  months  of  the  nonlinear  spin  down 
of  an  unstable  jet.  Weak  westward  recirculations  develop  in  the  lower  layer.  During 
the  next  28  months  of  integration,  the  jet  sharpens  again,  increases  its  amplitude  and 
develops  westward  recirculations  in  both  layers. 

A  set  of  linear  stability  analysis  of  broader  initial  basic  flows  indicates  the 
tendency  for  stronger  recirculations  when  the  basic  flow  is  strongly  baroclinic.  Figure 
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2-13  shows  the  results  of  eddy-mean  flow  interaction  for  the  jet  of  160  km  width.  No¬ 
tice  that  eddy  potential  vorticity  fluxes  act  to  decelerate  the  flow  in  the  upper  layer 
and  to  accelerate  the  flow  in  the  lower  layer.  The  residual  meridional  circulation  is 
responsible  for  diffusing  the  jet  in  the  upper  layer  and  for  inducing  the  westward 
retardation  in  the  lower  layer.  In  contrast  with  a  thinner  (100  km  wide)  jet,  the  am¬ 
plitudes  of  mean  zonal  accelerations  are  of  the  same  order  of  magnitude  in  both  layers. 
In  strongly  barotropic  flows,  the  tendency  for  westward  recirculations  is  absent.  As 
we  increase  the  width  of  the  jet,  the  flow  stabilizes. 

To  see  if  the  linear  mean-flow  interaction  alone  could  be  responsible  for  the  for¬ 
mation  of  the  recirculations  we  used  an  eddy-mean  flow  interaction  model,  developed 
by  Flierl  (personal  communication).  The  model  equations  are: 


=  (-I)'*1* +»;«;, 

(2.24) 

d£+tiM  +  vM  -o 

at  '  ’  dy 

(2.25) 

q2 

( 2  -  F\-  F2)x  =  V \qx  -  v2q2, 

(2.26) 

where  u  and  q  are  mean  zonal  geostrophic  velocity  vector  and  mean  potential  vortic- 
ity,  (u  ,  v  )T  and  q  are  perturbation  geostrophic  velocity  and  perturbation  potential 
vorticity,  x  is  residual  meridional  circulation  and  i  =  1, 2  is  layer  index. 

The  model  was  applied  for  a  basic  flow  of  the  form  u^y)  =  A,cosh2(y)  and 
converged  only  for  strongly  baroclinic  initial  basic  flows;  the  procedure  was  divergent 
for  strongly  barotropic  flows.  The  resulting  basic  flow  had  well  developed  westward 
flows.  Thus,  eddy-mean  flow  interactions  were  sufficient  to  explain  the  development 
of  recirculations  in  strongly  baroclinic  cases. 

This  approach  does  not  work  for  barotropic  flows,  since  strongly  nonlinear 
regimes  allowing  for  both  eddy-mean  flow  and  eddy-eddy  interactions  are  required 
for  recirculation  development. 
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Figure  2-14:  Eddy  -  mean  flow  partition  of  total  energy  vs  dimensional  time.  The 
dark-shaded  area  indicates  the  mean  flow  energy  part,  the  white  area  indicates  the 
eddy  energy  part.  Parameters:  ux  =  1,  u2  =  0,  (3  =  4.15,  7  =  10.18,  5  =  0.25,  U  = 
0.8  m  sec-1. 

2.3  Nonlinear  Analysis  of  a  Two-Layer  Model 

2.3.1  The  Flow  Evolution 

In  all  the  numerical  experiments,  flow  evolution  goes  through  two  separate  stages. 
During  the  first  stage,  the  eddy  energy  grows  linearly  at  the  expense  of  the  mean 
flow  energy  (Figure  2-14).  In  the  streamfunction  field  (Figure  2-15),  this  stage  corre¬ 
sponds  to  the  development  of  instabilities,  to  the  growth  of  meanders  and,  in  strongly 
barotropic  cases,  to  ring  formation.  The  second  stage  begins  when  the  eddy  energy 
reaches  its  maximum.  Eddies,  radiated  away  from  the  jet,  start  releasing  energy  back 
to  the  mean  flow  thus  restructuring  it.  The  equilibrated  flow  has  a  wider  and  weaker 
eastward  jet  surrounded  by  westward  flows  to  the  north  and  south. 

Total  integration  time  of  the  model  is  equivalent  to  30  months.  Total  net 
transport  is  constant  over  time  (Figure  2-16).  Both  integrated  eastward  and  westward 
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Figure  2-16:  Dimensional  total  net  (solid),  eastward  (dash-dot)  and  westward  (dash) 
transports  as  functions  of  dimensional  time.  Parameters:  ui  =  1,  u2  =  0,  ft  = 
4.15,  7  =  10.18,  6  =  0.25,  U  =  0.8  m  sec”1. 
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Figure  2-17:  Total  (solid),  potential  (dash),  upper  layer  kinetic  (dash-dot)  and  lower 
layer  kinetic  (dot)  energy.  Parameters:  U\  =  1,  u2  =  0,  (5  —  4.15,  7  =  10.18,  5  = 
0.25,  U  =  0.8  m  sec-1. 
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Figure  2-18:  Total  enstrophy  vs  dimensional  time.  Parameters:  U\  —  1,  u2  =  0,  /3  = 
4-15,  7  =  10.18,  6  =  0.25,  l/  =  0.8m  sec-1. 

transports  increase  linearly  for  the  first  7  months,  saturating  completely  after  15 
months. 

Since  the  model  has  two  sponge  layers  and  a  weak  numerical  filter,  total  energy 
slowly  decreases  with  time  (Figure  2-17).  The  maximum  decay  corresponds  to  the 
first  5  months  of  model  integration,  when  intensely  growing  eddies  radiate  away  from 
the  jet  and  enter  the  sponge  layers.  Once  eddies  start  releasing  energy  back  to  the 
mean  flow,  the  energy  decay  diminishes  substantially. 

Time  fluctuations  of  total  momentum  are  less  than  0.005%  of  the  total  initial 
field.  Total  enstrophy  decreases  with  time  by  about  5%  (Figure  2-18),  with  maximum 
decay  during  the  stage  of  eddy  energy  growth. 

2.3.2  The  Equilibrated  Jet  Structure 

This  sections  discusses  the  differences  in  the  equilibrated  jet  structures  for  baroclinic 
and  barotropic  initial  basic  flows.  The  initial  basic  velocity  is  chosen  such  that  ui  =  1 
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Figure  2-19:  Upper  (solid)  and  lower  (dash)  layer  mean  zonal  velocity,  Ui(y)  = 
2~  S—7T  Ui(x,  y)dx,  i  =  1, 2,  at  t=0  (a),  28  days  (b),  47  days  (c),  66  days  (d),  3  months 
(e),  4  months  (f),  6  months  (g),  30  months  (h).  Parameters:  U\  =  1,  u2  =  0,  /3  = 
4.15,  7  =  10.18,  8  =  0.25,  U  =  0.8  m  sec-1. 
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and  u2  varies  between  zero  and  0.99.  In  all  the  experiments,  the  initially  unstable  jet 
stabilizes,  and  its  structure  changes  drastically  with  time. 

Figure  2-19  shows  the  time  evolution  of  the  mean  zonal  velocity  for  the  case 
when  initially  there  is  no  flow  in  the  lower  layer  (a).  In  the  first  month  (b),  when 
eddies  grow  at  the  expense  of  the  mean  flow  energy,  the  upper  layer  jet  decreases  its 
amplitude  and  slightly  widens.  In  the  lower  layer,  a  weak  eastward  flow  develops  in 
the  middle  of  the  domain  surrounded  by  two  weak  westward  flows.  Between  1.5  (c) 
and  2  months  (d)  of  flow  evolution,  this  lower  layer  structure  strengthens,  while  the 
upper  layer  mean  zonal  velocity  loses  its  well  defined  jet  profile  due  to  the  still  growing 
eddy  field.  After  6  months  (g),  eddies  pump  energy  back  to  the  mean  flow,  and  the 
equilibrated  flow  (h)  consists  of  a  sharp  baroclinic  jet  (with  a  nonzero  amplitude  in 
the  lower  layer)  surrounded  by  barotropic  westward  flows  on  the  north  and  south. 
Figure  2-20  shows  the  mean  flow  evolution  for  strongly  barotropic  flows.  The  process 
is  similar  for  strongly  barotropic  basic  flows  (Figure  2-20) . 

As  the  nondimensional  lower  layer  velocity  amplitude,  u2,  varies  from  zero  to 
0.99,  the  final  maximum  of  the  mean  dimensional  zonal  velocity  decreases  from  0.88  to 
0.53  m  sec-1  in  the  upper  layer  and  increases  from  0.23  to  0.53  m  sec-1  in  the  lower 
layer  (Figure  2-21a).  Interestingly,  when  u2  =  0,  the  maximum  final  upper  layer 
velocity  is  larger  than  the  maximum  initial  upper  layer  velocity.  The  dimensional 
width  of  the  upper  layer  jet  increases  monotonically  from  180  to  220  km.  In  the  lower 
layer,  it  diminishes  from  140  to  120  km  as  u2  increases  from  0.0  to  0.5,  and  then  grows 
to  210  km  for  it2  >  0.5. 

2.3.3  Potential  Vorticity  Mixing 

The  equilibration  of  an  unstable  jet  can  be  viewed  through  the  process  of  potential 
vorticity  mixing  that  acts  to  remove  a  negative  potential  vorticity  gradients.  Such  a 
process  differs  for  baroclinic  versus  barotropic  jets.  In  the  following  discussion,  we 
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Figure  2-20:  Upper  layer  mean  zonal  velocity,  Ui(y )  =  ^  f*w  Ui(x,  y)dx,  at  t=0  (a), 
18  days  (b),  47  days  (c),  66  days  (d),  9  months  (e),  30  months  (f).  Parameters: 
iii  =  1,  u2  =  0.99,  j3  —  4.15,  7  -  10.18,  5  =  0.25,  U  =  0.8  m  sec-1. 
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Figure  2-21:  (a):  Maximum  of  dimensional  velocity,  U  uh  (  in  m  sec-1)  of  the  equi¬ 
librated  jet  in  the  upper  (solid)  and  lower  (dash)  layers  versus  the  nondimensional 
amplitude  of  the  initial  lower  layer  velocity,  u2.  Parameters:  iq  =  1,  /3  =  4.15,  7  = 
10.18,  8  =  0.25,  U  =  0.8  m  sec-1,  (b):  Dimensional  width  (in  km)  of  the  equilibrated 
jet  in  the  upper  (solid)  and  lower  (dash)  layer  versus  the  nondimensional  amplitude 
of  the  initial  lower  layer  velocity,  u2. 

contrast  two  examples:  {ui,u2)  =  (1,0)  which  corresponds  to  a  strongly  baroclinic 
jet,  and  («i,w2)  =  (1.0,0,99)  which  corresponds  to  a  strongly  barotropic  jet. 

For  u2  =  0,  both  the  upper  and  the  lower  layer  potential  vorticity  fields  have 
regions  of  negative  gradient  (Figure  2-22a).  After  equilibration,  the  lower  layer  poten¬ 
tial  vorticity  is  homogenized  in  a  wide  zone  under  the  axis  of  the  jet  (Figure  2-22b). 
Hence,  the  flow  is  only  marginally  stable  there.  Zero  g-gradient  indicates  that  there 
is  intense  mixing  in  the  lower  layer  in  the  center  of  the  domain.  Similar  homoge¬ 
nization  of  lower  layer  potential  vorticity  was  observed  in  a  two-layer  /3-plane  model 
of  Pedlosky  [28],  where  he  studied  weakly  nonlinear  equilibration  of  Phillips’  model 
with  weak  dissipation. 

In  the  rest  of  the  lower  layer  domain,  there  is  a  weak  positive  potential  vor¬ 
ticity  gradient  which  corresponds  to  planetary  /?.  In  the  upper  layer,  the  region  of 
the  eastward  jet  with  large  and  positive  g-gradient  is  surrounded  on  the  north  and 
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Figure  2-22:  Initial  (a)  and  final  (b)  mean  zonal  potential  vorticity,  ql(y)  = 
^  I-vQi(x,y)dx,  i  =  1,2,  for  the  upper  (solid)  and  lower  (dash)  layer.  Parameters: 
ux  =  1,  u2  =  0,  (3  =  4.15,  7  =  10.18,  8  =  0.25,  U  =  0.8  m  sec-1. 
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Figure  2-23:  Initial  (a)  and  final  (b)  mean  zonal  potential  vorticity  for  the  upper  layer, 
qi(y)  =  ^  f*7Tqi(x,y)dx.  Parameters:  u\  =  1,  u2  =  0.99,  /?  =  4.15,  7  =  10.18,  8  = 
0.25,  U  =  0.8  m  sec-1. 
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Figure  2-24:  (a):  Contours  of  the  upper  layer  potential  vorticity  at  t  =  47  days,  (b): 
Same  for  the  lower  layer.  Parameters:  u\  —  1,  u2  =  0,  (3  =  4.15,  7  =  10.18,  6  = 
0.25,  U  =  0.8  m  sec-1. 
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Figure  2-25:  (a):  Contours  of  the  upper  layer  potential  vorticity  at  t  =  66  days,  (b): 
Same  for  the  lower  layer.  Parameters:  U\  =  1,  =  0,  p  =  4.15,  7  =  10.18,  8  = 

0.25,  17  =  0.8  m  sec-1. 
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Figure  2-26:  (a):  Contours  of  the  upper  layer  potential  vorticity  at  t  =  4  months, 
(b):  Same  for  the  lower  layer.  Parameters:  ui  =  1,  U2  =  0,  f3  =  4.15,  7  =  10.18,  S  = 
0.25,  17  =  0.8  m  sec-1. 
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Figure  2-28:  Contours  of  the  upper  layer  potential  vorticity  at  t  =  18  days.  Parame¬ 
ters:  Mi  =  1,  u2  =  0.99,  /?  =  4.15,  7  =  10.18,  5  =  0.25,  U  —  0.8  m  sec-1. 
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Figure  2-29:  Contours  of  the  upper  layer  potential  vorticity  at  t  =  47  days.  Parame¬ 
ters:  ui  =  1,  u2  =  0.99,  /?  =  4.15,  7  =  10.18,  6  =  0.25,  U  =  0.8  m  sec-1. 
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Figure  2-30:  Contours  of  the  upper  layer  potential  vorticity  at  t  =  3  months.  Param¬ 
eters:  U\  =  1,  u2  =  0.99,  j3  =  4.15,  7  =  10.18,  6  =  0.25,  U  =  0.8  m  sec-1. 
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Figure  2-31:  Contours  of  the  upper  layer  potential  vorticity  at  t  =  30  months.  Pa¬ 
rameters:  Ui  =  1,  u2  —  0.99,  /?  =  4.15,  7  =  10.18,  6  =  0.25,  U  =  0.8  m  sec-1. 
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the  south  by  the  zones  of  westward  flows  with  homogenized  potential  vorticity.  Inter¬ 
estingly,  the  value  of  the  positive  q— gradient  decreased  only  by  10%  during  the  flow 
evolution,  which  indicates  very  weak  mixing  across  the  upper  layer  potential  vorticity 
front  and  intense  mixing  on  the  flanks  of  the  jet. 

For  the  case  of  a  strongly  barotropic  jet,  (ui,u2)  =  (1,0.99),  the  regions  of 
initial  negative  potential  vorticity  gradient  (Figure  2-23a)  are  removed  after  equili¬ 
bration  (Figure  2-23b).  The  final  q— structure  consists  of  a  zone  of  positive  potential 
vorticity  gradient  on  the  axis  of  the  jet  with  the  regions  of  homogenized  potential  vor¬ 
ticity  north  and  south  of  it.  In  contrast  with  a  strongly  baroclinic  case,  the  value  of 
positive  q— gradient  decreased  by  45%  ,  which  suggests  that  in  a  strongly  barotropic 
flow,  mixing  occurred  not  only  on  the  flanks  of  the  jet,  but  also  across  the  potential 
vorticity  front  on  the  axis  of  the  jet. 

The  differences  in  potential  vorticity  mixing  are  illustrated  in  Figures  2-24-2- 
31.  For  a  strongly  baroclinic  case,  (tq,  u2)  —  (1,  0),  a  large  positive  potential  vorticity 
gradient  in  the  upper  layer  jet  appears  as  an  efficient  barrier  to  mixing;  even  though 
meanders  grow  to  large  amplitudes,  the  front  stays  quite  coherent  and  breaks  only  in 
a  few  places  giving  birth  to  rings  (Figures  2-24,  2-25).  Therefore,  in  the  upper  layer, 
mixing  occurs  on  the  sides  of  the  jet;  potential  vorticity  filaments  get  strained  and 
stirred  on  each  side  of  a  meander.  Most  of  q— mixing  occurs  in  the  lower  layer,  where 
the  initial  potential  vorticity  gradient  is  negative.  After  4  months  of  integration  (Fig¬ 
ure  2-26),  the  amplitudes  of  the  upper  layer  meanders  decrease,  and  straining  of  the 
potential  vorticity  contours  in  the  lower  layer  continues.  Both  layers  are  dominated 
by  heton-like  structures.  The  equilibrated  field  (Figure  2-27)  is  characterized  by  ho¬ 
mogenized  lower  layer  potential  vorticity  under  the  axis  of  the  jet  and  by  a  strong 
q— front  in  the  upper  layer. 

In  the  strongly  barotropic  case  ( u2  =  0.99),  the  meanders  start  breaking  long 
before  they  propagate  across  the  domain  (Figure  2-28).  After  1  month  of  integration 
(Figure  2-29),  there  is  an  intense  potential  vorticity  mixing  across  the  front  and  on 
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its  sides.  The  front  becomes  completely  broken  after  3  months  (Figure  2-30).  Even 
though  initially  there  is  a  strong  potential  vorticity  gradient  in  the  middle  of  the  do¬ 
main,  which  should  serve  as  a  barrier  to  mixing,  the  eddy  field  is  so  intense  and  violent 
that  it  eventually  breaks  the  front.  This  contrasts  with  a  strongly  baroclinic  case  de¬ 
scribed  earlier.  The  equilibrated  potential  vorticity  field  (Figures  2-31)  has  a  positive 
gradient,  but  it  is  much  weaker  than  the  initial  one.  Two  regions  of  homogenized  q 
are  located  north  and  south  of  the  front. 

Energy  and  Momentum  Constraints 

During  a  free  spin  down,  why  does  the  flow  equilibrate  itself  by  developing  the  recircu¬ 
lations?  The  equilibration  process,  which  we  view  as  potential  vorticity  mixing,  must 
obey  the  conservation  of  both  momentum  and  energy.  One  should  remember  that 
although  the  model  is  unforced,  the  conservation  is  approximate  due  to  the  presence 
of  a  numerical  filter  and  the  sponge  layers. 


Figure  2-32:  A  schematic,  illustrating  a  thought  experiment.  Initial  velocity  profile 
(left)  and  equilibrated  velocity  profile  (right). 


To  illustrate  the  importance  of  the  above  constraints,  let  us  consider  the  fol¬ 
lowing  thought  experiment.  Suppose  that  an  initial  barotropically  unstable  “point” 
eastward  jet  has  a  velocity  amplitude  of  10  units  and  the  kinetic  energy  of  ^  units 
(Figure  2-32,  left).  Due  to  potential  vorticity  mixing,  the  equilibrated  flow  can  be 
represented  by  a  broadened  jet  with  velocity  amplitudes  {02,01,00,01,02},  where 
o2  <  Oi  <  a0  (Figure  2-32,  right).  Note  that  we  make  no  assumptions  about  the  sign 
of  the  amplitude  components,  a2.  The  conservation  of  momentum  and  energy  imply: 

Oo  +  2(ai  -f-  a2)  =  10, 

-  „  (2.27) 

Oq  +  2(a2  +  a2)  —  102. 

Figure  2-33  shows  the  solutions  of  the  above  equations,  when  the  center  amplitude, 
o0, is  considered  as  an  independent  parameter  varying  between  0  and  10.  This  simple 
experiment  shows  that  the  equilibrated  jet  structure  must  have  negative  velocity 
amplitudes,  a2,  on  the  flanks  of  the  equilibrated  barotropic  jet  in  order  to  conserve 
both  energy  and  momentum. 


Figure  2-33:  The  solutions  of  system  of  equations  (2.27)  conserving  momentum  and 
energy.  The  velocity  amplitudes  <n  (solid  line)  and  a2  (dashed  line)  are  shown  as 
functions  of  central  velocity  amplitude,  ao- 
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Figure  2-34:  Zonally  averaged  relative  vorticity,  ^\yy,  (solid),  stretching  term,  F(l)  * 
(d/2  —  ^i),  (dash)  and  planetary  potential  vorticity,  f3y,  (dash-dot)  in  the  upper  layer 
at  the  initial  (a)  and  final  (b)  time.  Similar  terms  for  the  lower  layer  at  the  initial  (c) 
and  final  (d)  time.  Parameters:  U\  =  1,  u2  =  0,  (3  =  4.15,  7  =  10.18,  <5  =  0.25,  U  — 
0.8  m  sec-1. 


67 


Figure  2-35:  Zonally  averaged  relative  vorticity,  lw,  (solid),  stretching  term,  F{  1)  * 
(^2-^1),  (dash)  and  planetary  potential  vorticity,  /fy,  (dash-dot)  in  the  upper  layer 
at  the  initial  (a)  and  final  (b)  time.  Parameters:  ux  =  1,  u2  =  0.99,  /3  =  4.15,  7  = 
10.18,  S  =  0.25,  U  —  0.8  m  sec-1. 


Figure  2-36:  Final  ip  -  q  dependence  for  the  upper  layer  (a)  and  lower  layer  (b) 
when  U\  =  1,  u-i  =  0,  f3  =  4.15,  7  =  10.18,  5  =  0.25,  U  =  0.8  m  sec-1. 

Potential  Vorticity  Components 

The  changes  in  the  individual  components  of  potential  vorticity  fields  are  presented 
in  Figures  2-34  and  2-35.  In  strongly  baroclinic  cases,  although  the  stretching  term, 
(— 1  yFi(ipi  —  ip2),  i  =  1,2,  varies  somewhat,  most  of  the  changes  come  from  the 
zonal  relative  vorticity,  ipiyy,  in  the  center  of  the  domain.  Therefore,  changes  in 
the  meridional  structure  of  the  flow  are  mostly  responsible  for  removing  the  negative 
values  of  potential  vorticity  gradient.  In  the  recirculation  regions,  the  relative  vorticity 
term  is  zero  and  the  stretching  term  is  constant  in  y. 

ip  —  q  Relationship  for  Equilibrated  Jet 

Figure  2-36a  shows  three  different  regions  of  final  ip  —  q  relationship  in  the  upper  layer 
for  the  case  of  strongly  baroclinic  initial  basic  flow.  The  axis  of  the  jet  corresponds 
to  the  linear  and  negative  ip  —  q  dependence.  In  the  zones  of  recirculations,  where 
the  potential  vorticity  gradient  is  zero,  the  ip  —  q  function  is  horizontal.  The  third 
zone,  where  ip  —  q  is  positive  and  linear,  is  in  the  area  to  the  north  and  south  of  the 
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Figure  2-37:  Final  ip  -  q  dependence  for  the  upper  layer  (a)  and  lower  layer  (b) 
when  =  1,  «2  =  0.99,  0  =  4.15,  7  =  10.18,  8  =  0.25,  U  =  0.8  m  sec”1. 

westward  flows.  Here,  the  potential  vorticity  gradient  is  equal  to  0.  In  the  lower  layer 
(Figure  2-36b),  both  in  the  region  of  eastward  jet  and  westward  recirculations,  the 
ip  -q  dependence  is  horizontal.  For  strongly  barotropic  initial  basic  flow  (Figure  2- 
37),  the  final  ip-q  relationship  is  similar  to  that  of  the  upper  layer  strongly  baroclinic 
case. 

2.3.4  Relative  Strength  of  the  Recirculations 

Dependence  on  Baroclinic  Structure 

To  understand  how  the  strength  of  the  recirculations  depends  on  the  baroclinic  struc¬ 
ture  of  the  jet  for  a  given  /?,  we  conducted  a  series  of  experiments  in  which  the 
amplitude  of  the  nondimensional  upper  layer  basic  velocity  was  fixed  at  1  and  the 
amplitude  of  the  nondimensional  lower  layer  basic  velocity  was  varied  between  0  and 
0.99. 

Figure  2-38  shows  total  dimensional  net,  eastward  and  westward  transports  for 
the  equilibrated  jet  as  functions  of  the  lower  layer  nondimensional  velocity  amplitude, 
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Figure  2-38:  Total  net  (o),  eastward  (*)  and  westward  (□)  transports  (in  Sv)  for  the 
equilibrated  jet  as  functions  of  the  amplitude  of  the  nondimensional  lower  layer  basic 
velocity,  u2.  The  upper  layer  nondimensional  velocity  amplitude,  iq,  was  fixed  at  1.0. 
Parameters:  (3  =  4.15,  7  =  10.18,  S  =  0.25,  U  =  0.8  m  sec-1.  Dashed  lines  show  the 
observed  transport  in  the  Gulf  Stream  (  150  Sv  )  and  in  the  recirculations  ( -80  Sv  ). 

u2,  for  P  =  4.15.  Both  eastward  and  westward  transports  were  not  monotonic.  For 
u2  <  0.3,  the  eastward  transport  decreased  from  250  Sv  to  200  Sv,  whereas  the 
westward  transport  diminished  from  180  Sv  to  30  Sv.  For  u2  >  0.3,  they  both 
increased  to  500  Sv  and  110  Sv  correspondingly.  Note  that  the  observed  transports 
in  the  Gulf  Stream  and  its  recirculations  correspond  to  150  Sv  and  80  Sv  ( Hogg  [12]). 

Similar  non-monotonic  dependence  of  the  recirculations  strength  on  the  lower 
layer  nondimensional  velocity  amplitude  was  observed  for  three  different  values  of 
nondimensional  ft  (Figure  2-39)  corresponding  to  three  different  velocity  scales:  U  = 
0.6, 0.8  and  1.0  m  sec-1.  The  minimum  in  the  strength  of  the  recirculations  was 
attained  at  u2  =  0.3  for  P  =  4.15  and  (3  =  5.53  and  at  u2  =  0.5  for  P  =  3.32. 
Therefore,  the  recirculations  were  strongest  for  strongly  baroclinic  initial  basic  flows, 
weakest  for  flows  with  intermediate  baroclinicity  and  of  medium  strength  for  strongly 
barotropic  flows. 
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Figure  2-39:  The  relative  strength  of  the  recirculations,  defined  by  the  ratio  of  total 
westward  to  eastward  transports,  as  a  function  of  the  lower  layer  nondimensional 
velocity  amplitude,  it2.  Parameters:  ux  =  1,  7  =  10.18,  5  =  0.25.  The  dashed,  solid 
and  dashed-dotted  lines  correspond  to  (3  =  3.32,  U  —  1  m  sec-1,  /5  =  4.15,  U  =  0.8 
m  sec  1  and  /?  =  5.53,  U  ~  0.6  m  sec-1  respectively. 


Figure  2-40.  (a).  The  ratio  of  the  total  westward  to  the  total  eastward  transport 
as  a  function  of  a  baroclinic  parameter,  .  (b):  The  dimensional  total  net 

(solid  line),  westward  (dashed-dotted  line)  and  eastward  transports  (dashed  line)  as 
functions  of  a  baroclinic  parameter.  Parameters:  (3  =  4.15,  7  =  10.18,  5  =  0.25. 
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Figure  2-41:  Maximum  linear  growth  rate  as  a  function  of  a  baroclinic  parameter, 
f1uI+'f2u1  •  Parameters:  f3  =  4.15,  72  =  103.75,  8  =  0.25. 


However,  the  above  experiments  differed  in  their  total  energy  levels.  To  see  if 
the  existence  of  the  local  minimum  in  the  recirculation  strength  could  be  related  to 
changes  in  the  total  energy  level  or  in  the  partition  of  kinetic  to  potential  energies, 
a  set  of  experiments  was  conducted  where  the  total  energy  was  fixed  and  only  the 
baroclinicity  of  the  basic  flow  was  changed.  The  latter  was  defined  as  the  ratio  of 
the  basic  baroclinic  velocity,  u\  -  u2,  to  the  basic  barotropic  velocity,  Fxu2  +  F2u\. 
Figure  2-40a  shows  the  resulting  relative  strength  of  the  recirculations.  There  is 
still  a  minimum  reached  at  the  value  of  the  baroclinic  parameter  that  corresponds  to 
«i  =  1.0  and  u2  =  0.3.  The  corresponding  partition  of  kinetic  to  potential  energies  was 
a  monotonic  function  of  the  baroclinic  parameter.  The  total  dimensional  transport 
(Figure  2-40b)  also  indicates  that  minimum  absolute  value  of  the  westward  transport 
is  attained  at  Mi  =  1.0  and  u2  —  0.3. 

Therefore,  a  non-monotonic  behavior  of  the  relative  strength  of  the  recircula¬ 
tions  is  not  related  to  the  total  energy  level  or  to  the  kinetic-potential  energy  partition. 
Instead,  it  is  determined  solely  by  the  baroclinic  structure  of  the  initial  basic  flow. 
Linear  stability  analysis  reported  in  section  2.2.3  showed  that  the  growth  rate  was  a 
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non-monotonic  function  of  the  lower  layer  nondimensional  velocity,  u2,  for  a  fixed  Ui 
and  (5  (Figure  2-6).  The  shape  of  the  growth  rate  curves  indicated  the  importance 
of  a  baroclinic  mode  when  u2  was  close  to  zero  and  of  a  barotropic  mode  when  u2 
was  close  to  uv  For  the  intermediate  values  of  u2,  the  growth  rate  was  minimum,  and 
both  modes  had  similar  growths.  Also,  in  the  experiments  with  the  fixed  total  energy 
level,  maximum  growth  rate  was  a  non-monotonic  function  of  a  baroclinic  parameter, 
(Figure  2-41). 

Linear  stability  analysis  also  showed  that  recirculation  development  due  to 
linear  wave-mean  flow  interaction  was  possible  only  for  strongly  baroclinic  flows.  An 
eddy-mean  flow  interaction  model  applied  to  a  strongly  baroclinic  basic  flow  produced 
recirculations  in  both  layers.  The  same  procedure  applied  to  a  strongly  barotropic 
flow  did  not  result  in  flow  equilibration.  Therefore,  for  strongly  baroclinic  basic 
flows,  eddy-mean  flow  interaction  alone  can  generate  recirculations  in  both  layers 
through  residual  meridional  circulation.  Thus,  the  results  of  the  nonlinear  calculations 
showing  a  decrease  in  the  relative  strength  of  the  recirculations  with  an  increase  in 
u2  for  strongly  baroclinic  jets  are  in  agreement  with  the  results  of  the  linear  stability 
analysis. 

On  the  other  hand,  recirculations  developed  in  the  process  of  a  nonlinear  equi¬ 
libration  of  a  strongly  barotropic  flow.  Clearly,  nonlinear  eddy-eddy  and  eddy-mean 
flow  interactions  were  necessary  for  their  generation.  Therefore,  the  principle  mech¬ 
anism  leading  the  recirculation  development  is  different  in  the  case  of  strongly  baro¬ 
clinic  and  strongly  barotropic  flows. 

In  the  present  study,  we  did  not  perform  a  weakly-nonlinear  analysis  of  the 
basic  state  flows.  However,  Flierl  and  Meacham  (personal  communication)  showed 
that  in  a  contour  dynamical  model  of  a  two-layer  cusp  jet  with  three  potential  vorticity 
fronts  in  each  layer,  the  coefficient  in  the  amplitude  equation  changes  sign  for  the 
intermediate  values  of  the  baroclinic  parameter,  indicating  that  the  flow  does  not 
saturate.  Therefore,  equilibration  requires  a  strongly  nonlinear  regime.  Similarly,  in 
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Figure  2-42:  The  ratio  of  the  upper  to  the  lower  layer  westward  transports  for  the 
equilibrated  jet,  as  a  function  of  the  lower  layer  nondimensional  velocity  amplitude, 
u2.  Parameters:  iq  =  1,7  =  10.18,  6  =  0.25.  The  dashed,  solid  and  dashed-dotted 
lines  correspond  to  /3  =  3.32,  V  —  1  m  sec-1,  /3  =  4.15,  U  =  0.8  m  sec-1  and 
P  =  5.53,  U  =  0.6  m  sec-1  respectively. 

a  case  of  continuous  potential  vorticity  gradient,  a  weakly-nonlinear  analysis  applied 
to  basic  state  flows  with  intermediate  baroclinicity  might  not  necessarily  give  any 
insight  into  the  equilibration  process. 

Figure  2-42  presents  the  ratio  of  the  upper  to  the  lower  layer  westward  trans¬ 
ports  of  the  equilibrated  jet,  as  a  function  of  the  lower  layer  nondimensional  velocity 
amplitude,  u2,  for  three  different  values  of  nondimensional  /3.  The  ratio  varied  be¬ 
tween  0.15  and  0.24  with  maximum  values  reached  at  the  intermediate  values  of  u2 
and  minimum  values  attained  when  u2  =  0.  Recalling  that  the  layer  depths  were  cho¬ 
sen  as  1  km  and  4  km  suggests  that  lower  layer  recirculations  were  always  stronger 
than  upper  layer  recirculations. 


Dependence  on  /3 

The  relative  strength  of  the  recirculations,  defined  as  the  ratio  of  the  total  westward 
to  the  total  eastward  transport  of  the  equilibrated  jet,  increased  with  a  decrease  in 
the  value  of  a  nondimensional  /3-parameter  (Figure  2-39).  This  result  is  in  agreement 
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Figure  2-43:  Schematic  showing  mean  potential  vorticity  for  two  values  of  nondimen- 
sional  p  :  Pi  >  p2.  Light  shaded  area  indicates  mixed  q  for  p2,  dark  shaded  area 
indicates  mixed  q  for  Pi. 

with  the  findings  by  Jayne  et  al.  [18],  who  studied  the  statistically  steady  state  of 
the  eastward  jet  and  its  recirculations  driven  by  the  unstable  jet  inflow  on  the  eastern 
boundary  and  by  the  stable  jet  outflow  though  the  western  boundary. 

The  increase  in  the  recirculation  strength  with  a  decrease  in  nondimensional  p 
can  be  easily  interpreted  by  viewing  the  equilibration  process  as  the  one  of  potential 
vorticity  mixing.  Figure  2-43  shows  schematic  potential  vorticity  profiles  for  two 
different  values  of  nondimensional  P  \  P\  >  p2.  After  equilibration,  the  regions  of 
negative  potential  vorticity  gradient  are  removed.  Recirculations  appearing  north 
and  south  of  the  jet  correspond  to  the  regions  with  homogenized  potential  vorticity. 
Let  us  assume  here  that  total  enstrophy  is  conserved  and  that  the  potential  vorticity 
gradient  on  the  axis  of  the  jet  does  not  change  during  equilibration.  Then,  the 
following  areas  must  be  equal  (see  Figure):  A\  =  A\  and  A\  =  A\.  Therefore,  the  size 
of  the  homogenized  region  is  smaller  for  P  =  /A,  which  indicates  smaller  recirculations. 
Hence,  conservation  of  enstrophy  and  symmetry  of  the  flow  imply  that  the  region 
of  homogenized  potential  vorticity  for  /fi  will  be  smaller  than  that  for  p2.  In  the 
present  model,  enstrophy  is  not  exactly  conserved  due  to  the  presence  of  weak  sponge 
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layers  and  a  numerical  filter;  it  decreases  by  less  than  5%  during  spin  down.  Also, 
the  potential  vorticity  gradient  corresponding  to  the  jet  axis  also  changes  depending 
on  the  baroclinicity  of  the  basic  flow.  However,  the  above  argument  illustrates  the 
importance  of  potential  vorticity  mixing  and  approximate  enstrophy  conservation  as 
constraints  on  the  recirculation  size. 


2.4  Results  for  the  Reduced- Gravity  Model 

To  understand  the  impact  of  the  barotropic  instability  on  the  mean  flow  and  its 
possible  role  in  generating  the  recirculations,  a  linear  stability  analysis  was  performed 
in  a  framework  of  a  reduced-gravity  model.  The  latter  can  be  considered  as  a  limiting 
case  of  a  two  layer  model,  described  in  section  2.1,  when  the  depth  of  the  lower  layer 
is  infinite,  i.e.  S  =  DxjD2  =  0.  The  upper  layer  is  dynamically  active  and  is  governed 
by  the  following  equation: 

^  =0,  (2.28) 

where  ’J'(y)  =  -  Jy  U(y')dy'  and  Q(y )  =  fiy  -  72\k  +  are  the  basic  flow  stream- 
function  and  potential  vorticity  fields,  </>  =  (j){x ,  y ,  t )  and  q  —  —y2<\>  +  are  the  per¬ 
turbation  streamfunction  and  potential  vorticity  fields.  The  model  is  defined  by  the 
following  nondimensional  parameters:  the  inverse  deformation  radius,  7  =  L/Ldef, 
and  nondimensional  beta-parameter,  /3  =  /3dimL2/U.  Basic  state  velocity  is  given  by 
Equation  (2.5). 

Although  the  structure  of  the  basic  mode  is  baroclinic,  only  the  mechanism  of 
barotropic  instability  is  present  in  the  reduced-gravity  model.  Again,  we  look  for  a 
zonally  periodic  solution  of  the  form 

<f>(x,  y,  t)  =  Real($(y)  ett<x-rf>),  (2.29) 
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figure  2-44:  The  maximum  growth  rate,  kcu  (solid  line)  and  the  corresponding  real 
part  of  the  phase  speed,  cr,  (denoted  by  circles)  as  functions  of  a  horizontal  wave 
number,  k,  in  the  reduced-gravity  model  for  7  =  9.11,  Ldef  =  44  km  (a),  7  = 
6.27,  Ldef  =  65  km  (b),  7  =  4.07,  Ldef  =  100  km  (c)  and  7  =  0,  Ldef  =  00  (d). 
Parameters:  (5  =  4.15,  w  =  8.14. 


Figure  2-45:  Zonal  mean  acceleration,  ||,  (a),  northward  potential  vorticity  flux,  vq, 
(b)  and  residual  meridional  circulation,  v*,  (c).  Parameters  of  the  reduced-gravity 
model:  Ui  =  1,  (3  =  4.15,  7  =  9.11,  w  =  8.14. 
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Figure  2-46:  Zonal  mean  acceleration,  in  a  reduced-gravity  model.  Parameters: 
Ui  =  1,  P  —  4.15,  7  =  0,  w  =  8.14. 

where  &  is  a  horizontal  wave  number,  and  c  is  the  phase  speed  of  the  disturbance. 
Then  the  linearized  equation  of  motion  is  given  by 

(C/(y)-c)(~^-A:2$-72$)  +  $^  =  0,  (2.30) 

subject  to  no-normal-flow  boundary  conditions 

$  (y  =  tt)  =  $  (y  =  -tt)  =  0.  (2.31) 

As  in  a  two  layer  model,  the  above  problem  is  reduced  to  the  generalized  eigenvalue 
problem: 


A$  =  c£$,  (2.32) 

where  matrices  A  and  B  are  given  by  (2.12)  and  (2.13). 

Figure  2-44  shows  the  maximum  growth  rates  and  the  real  parts  of  the  phase 
speed  for  four  values  of  nondimensional  7,  which  correspond  to  the  following  values 
of  the  deformation  radius:  44  km,  65  km,  100  km  and  infinite.  There  is  substantial 
difference  in  the  maximum  growth  rates.  When  the  deformation  radius  is  44  km, 
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the  flow  is  only  weakly  unstable,  and  there  is  long  wave  cutoff.  As  the  deformation 
radius  decreases,  the  long  waves  become  weakly  unstable.  The  most  unstable  wave 
increases  from  256  to  285  km  with  a  decrease  in  7.  The  corresponding  real  parts  of 
the  phase  speed  are  similar  for  various  7.  As  the  deformation  radius  increases,  short 
waves  become  more  unstable. 

The  TEM  analysis  shows  the  diffusion  and  the  rectification  of  the  flow  in  the 
middle  of  the  domain  and  the  development  of  westward  flows  on  the  flanks  of  the 
jet  when  Ldef  =  44  km  (Figure  2-45 (a)).  The  diffusion  of  the  jet  is  due  to  the  eddy 
potential  vorticity  fluxes  (b),  whereas  westward  rectification  is  due  to  the  residual 
meridional  circulation  (c).  In  the  case  of  an  infinite  deformation  radius,  Ldef  —  00, 
the  reduced-gravity  model  becomes  purely  barotropic.  Similar  to  a  two  layer  case,  the 
TEM  analysis  does  not  show  any  tendency  for  the  formation  of  westward  rectifications 
(Figure  2-46);  the  eddy  potential  vorticity  fluxes  have  a  purely  diffusive  effect  on  the 
mean  flow. 

The  numerical  experiments  with  the  fully  nonlinear  reduced-gravity  model 
show  that  the  relative  strength  of  the  recirculations,  defined  by  the  ratio  of  the  total 
westward  to  the  total  eastward  transports,  is  a  monotonically  decreasing  function 
of  a  nondimensional  inverse  deformation  radius,  7,  (Figure  2-47).  When  7  =  0, 
which  corresponds  to  an  infinite  dimensional  deformation  radius,  the  total  westward 
transport  is  20%  of  the  total  eastward  transport,  which  compares  well  with  the  results 
of  the  two-layer  model  for  a  strongly  barotropic  basic  flow.  When  the  deformation 
radius  was  chosen  to  be  44  km  (7  =  9.11),  the  strength  of  the  recirculations  was 
negligible. 

Therefore,  in  the  case  of  the  reduced-gravity  model,  the  tendency  for  the  for¬ 
mation  of  the  recirculations  during  the  stage  of  linear  growth  is  misleading  in  that  it 
predicts  the  development  of  recirculations  for  large  7  and  absence  of  the  recirculations 
for  small  7.  The  initial  basic  flow  is  only  slightly  supercritical  for  large  7,  and  therefore 
there  are  no  strongly  nonlinear  wave-mean  and  wave-wave  interactions  which  would 
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Figure  2-47:  The  ratio  of  westward  to  eastward  transport  in  the  equilibrated  jet  vs 
nondimensional  7  in  a  reduced-gravity  model.  Parameters:  u\  =  1  (3  =  A  15  w  — 
8.14.  ’  . 

be  necessary  for  the  development  of  the  recirculations  by  barotropic  instability.  In 
the  case  of  large  7,  the  basic  flow  is  strongly  supercritical,  and  hence  the  nonlinear 
eddy  field  is  strong  enough  to  produce  westward  rectifications. 

2.5  Conclusions 

The  formation  of  westward  recirculations  is  studied  in  the  model  of  a  free  nonlinear 
spin  down  of  zonally  symmetric,  quasi-geostrophic,  unstable  jets.  The  model  ( Flierl , 
personal  communication)  has  weak  sponge  layers  and  a  numerical  filter,  which  imply 
only  weak  dissipation  of  energy  and  enstrophy. 

The  recirculations  are  robust  features  of  a  two-layer  model,  which  includes 
the  mechanisms  of  both  baroclinic  and  barotropic  instabilities.  The  relative  strength 
of  the  recirculations,  defined  as  the  ratio  of  the  total  westward  to  the  total  eastward 
transport,  is  a  function  of  the  nondimensional  parameters  of  the  model.  It  decreases  as 
the  value  of  nondimensional  /?  increases.  This  decrease  in  the  recirculation  strengths 
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is  due  to  a  decrease  in  the  zone  of  homogenized  potential  vorticity,  which  in  turn  is 
due  to  (approximate)  enstrophy  conservation. 

When  P,  7  and  total  energy  level  are  fixed,  the  strength  of  the  recirculations  is 
maximum  for  strongly  baroclinic  flows  and  minimum  for  intermediate  values  of  a  baro- 
clinic  parameter,  where  ui  is  i- th  layer  nondimensional  velocity  amplitude 

and  Fi  is  i — th  layer  Froude  number.  In  strongly  baroclinic  cases,  eddies  marginally 
equilibrate  the  flow  and  create  a  wide  zone  of  homogenized  potential  vorticity  in  the 
lower  layer.  This  result  is  in  agreement  with  the  study  of  a  weakly-nonlinear  equili¬ 
bration  of  Phillips’  model  with  weak  dissipation  ( Pedlosky  [28]).  In  the  upper  layer,  a 
positive  potential  vorticity  gradient  on  the  axis  of  the  jet  decreases  by  less  than  10%, 
thus  appearing  as  an  efficient  barrier  to  mixing.  Most  vigorous  mixing  occurs  on  the 
edges  of  the  jet  in  the  upper  layer  and  under  the  axis  of  the  jet  in  the  lower  layer. 

As  the  jet  becomes  more  barotropic,  a  positive  potential  vorticity  gradient 
develops  on  the  axis  of  the  jet  in  the  lower  layer,  reflecting  a  more  stable  flow  structure. 
In  strongly  barotropic  jets,  there  is  vigorous  mixing  both  across  the  positive  potential 
vorticity  gradient  and  on  the  edges  of  the  jet,  leading  to  a  decrease  on  a  values  of  the 
upper  layer  potential  vorticity  gradient  on  the  axis  of  the  jet  by  more  than  45%. 

The  linear  stability  analysis  applied  to  a  two-layer  jets  shows  that  in  strongly 
baroclinic  cases,  the  lower  layer  recirculations  are  driven  by  the  eddy  heat  fluxes 
through  the  residual  meridional  circulation.  In  the  upper  layer,  the  divergence  of  eddy 
potential  vorticity  fluxes  tends  to  simply  “diffuse”  the  jet  during  the  stage  of  linear 
growth.  Nonlinear  eddy-eddy  and  eddy-mean  flow  interactions  are  further  important 
in  inducing  the  upper  layer  recirculations.  The  linear  analysis  shows  that  the  strength 
of  the  residual  meridional  circulations  increases  with  an  increase  in  the  width  of  basic 
flow.  That  is  why  a  “mean  flow”  model  ( Flierl ,  personal  communication)  successfully 
reproduces  recirculations  when  applied  to  strongly  baroclinic  basic  flows. 

For  strongly  barotropic  two-layer  basic  flows,  which  are  strongly  supercritical, 
the  linear  eddy  potential  vorticity  fluxes  tend  to  diffuse  the  jet  in  both  layers  without 
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any  tendency  for  the  formation  of  westward  flows  on  the  flanks.  As  the  basic  flow 
becomes  wider,  it  stabilizes.  Therefore,  strongly  nonlinear  regimes,  which  allow  for 
both  eddy-mean  flow  and  eddy-eddy  interactions,  are  required  for  the  emergence  of 
recirculations  for  strongly  barotropic  basic  flows. 

Thus,  although  the  recirculations  develop  in  a  two-layer  model  for  barotropic 
and  baroclinic  flows,  processes  leading  to  their  generation  are  different.  In  baro- 
clinic  case,  eddy-mean  flow  interaction  is  powerful  enough  to  drive  recirculations.  In 
barotropic  cases,  strongly  nonlinear  eddy-mean  flow  and  eddy-eddy  interactions  are 
required  for  the  recirculations  development.  The  existence  of  two  different  mecha¬ 
nisms  is  responsible  for  a  non-monotonic  dependence  of  the  recirculation  strength  on 
the  baroclinic  parameter  on  nonlinear  experiments. 

In  the  reduced-gravity  model,  which  has  just  the  mechanism  of  barotropic 
instability,  recirculations  form  as  a  result  of  a  nonlinear  spin  down  of  an  unstable  basic 
flow  only  for  the  case  of  an  infinite  dimensional  deformation  radius.  For  deformation 
radii  smaller  than  65  km,  the  total  westward  transport  is  less  than  5%  of  the  total 
eastward  transport. 

The  results  of  the  linear  stability  analysis  are  misleading  since  they  imply  a 
weak  tendency  for  the  formation  of  the  recirculations  for  small  dimensional  deforma¬ 
tion  radii  and  the  absence  of  the  recirculations  for  an  infinite  dimensional  deformation 
radius.  The  linear  stability  analysis  is  not  helpful  since  the  basic  flow  is  only  slightly 
supercritical  for  large  7,  implying  a  weak  eddy  field. 
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Chapter  3 


Colliding  Jets  Model  Formulation 

3.1  Introduction 

In  the  following  Chapters,  we  investigate  the  problem  of  low  frequency  variability  and 
associated  multiple  dynamical  regimes  resulting  from  the  internal  dynamics  of  an  un¬ 
stable  quasi-geostrophic  jet  and  its  recirculations.  Particularly,  we  raise  the  following 
questions:  Are  low  frequency  variability  and  multiple  dynamical  regimes  an  inherent 
part  of  the  internal  dynamics  of  a  quasi-geostrophic  jet  and  its  recirculations?  Are 
there  any  differences  between  barotropic  and  baroclinic  instabilities  as  mechanisms 
affecting  low  frequency  variability?  Which  dynamical  regimes  are  associated  with 
low  frequency  variability  and  which  precesses  are  responsible  for  transitions  between 
them?  Are  there  differences  between  low  frequency  variability  and  multiple  dynami¬ 
cal  regimes  in  the  globally  forced  (wind-driven)  models  and  boundary-forced  models, 
which  we  use  in  a  present  study? 

To  study  the  above  questions,  we  designed  regional  reduced-gravity  and  two 
layer  quasi-geostrophic  models  that  are  forced  by  prescribed  inflow  and  outflow  bound¬ 
ary  conditions;  dissipation  is  the  only  active  forcing  term  in  the  equations  of  motion. 
The  open  boundaries  represent  the  border  between  the  “outside”  region  of  slow  Sver¬ 
drup  interior,  where  the  wind  forcing  plays  a  crucial  role,  and  the  “interior”  region 
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of  the  western  boundary  current,  its  eastward  extension  and  its  recirculations,  where 
the  internal  dynamics  of  the  jet  is  dominant.  In  our  model,  the  role  of  wind-forcing  is 
limited  to  defining  a  particular  distribution  of  potential  vorticity  and  streamfunction 
along  the  inflow  and  outflow  boundaries,  a  distribution  that  is  fixed  in  time.  Notice 
that,  in  the  wind  forced  models  which  exhibit  multiple  dynamical  regimes,  the  in¬ 
flow  into  the  jet  can  change,  since  Sverdrup  balance  is  not  exact;  such  changes  may 
influence  the  dynamical  behavior. 


3.2  Model  Formulation 

3.2.1  Governing  Equations 

The  numerical  model  is  a  discretized  version  of  the  time-dependent  two-layer  quasi- 
geostrophic  potential  vorticity  equations  on  a  ft  -  plane: 

+  JtyuQi)  =  Vi  +  Si  in  ft,  (3.1) 

subject  to  the  boundary  conditions: 

A  =  ^i,  V2A  =  A,  VVi  =  Mi  in  ft.  (3.2) 

Here,  A  =  ipfix,  y )  is  a  quasi-geostrophic  streamfunction,  q{  =  qfix,  y )  =  I- 

—  is  potential  vorticity,  and  i  =  1,  2  is  a  layer  index.  The  rectangular  model 
domain  is  defined  by  ft  U  ft,  where  fi  is  the  set  of  all  the  internal  points  and  0  is  the 
set  of  all  the  boundary  points.  On  the  boundaries,  the  streamfunction,  its  Laplacian 
and  the  square  of  the  Laplacian  are  prescribed  as  known  functions  vfq,  £*•  and  Mj. 
The  dissipation  is  provided  by  biharmonic  diffusion  and  by  bottom  or  interfacial 
(in  the  case  of  the  reduced-gravity  formulation)  friction,  so  that  Vx  = 
and  V 2  =  -A2V2ip2  ~  AV6^2-  In  addition,  the  dissipation  St  in  the  sponge  layers, 
of  which  details  are  presented  later,  prevents  part  of  the  outgoing  radiation  from 
reflecting  back  onto  the  interior  domain. 
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3.2.2  Dimensional  Scales 


The  dimensional  scales  are  chosen  as  follows.  The  horizontal  length  scale  is  L  = 
2400  km,  the  advective  velocity  scale  is  U  =  1  m  sec-1,  which  implies  an  advective 
time  scale  T  =  L/U  =  27.8  days.  The  layer  thicknesses  in  the  reference  run  are 
Di  —  600  m  and  D2  —  4  km.  The  Coriolis  parameter  is  f0  =  7.25  •  10-5  sec-1,  its 
meridional  gradient  is  f3dim  =  1.97T0-11  m-1  sec-1.  The  baroclinic  deformation  radius 
is  chosen  as  Ldef  =  47  km.  Reduced  gravity  is  g  =  0.02  m  sec-2.  The  dimensional 
biharmonic  diffusion  coefficient,  ^4gim,  equals  either  8-1010  m4  sec-1  or  8-108  m4  sec-1. 
The  values  of  the  dimensional  bottom  friction  or  interfacial  friction,  Afm ,  vary. 

3.2.3  Model  Domain 

The  model  domain  is  a  rectangular  basin  defined  as 

OuH  =  {{x,y)  |  xw  <  x  <  xe,  ys  <  y  <  yn}  ,  (3.3) 

where  xe  =  —xw  =  Lx/2Ly,  yn  —  —ys  =  0.5.  The  dimensional  west-east  and  south- 
north  extents  of  the  domain  are  Lx  =  4667  km  and  Ly  =  2400  km.  The  model  has  a 
uniform  grid  with  257  by  129  grid  points,  which  implies  a  horizontal  resolution  of  18.5 
km.  High  resolution  runs  with  513  by  257  grid  points  and  a  corresponding  horizontal 
resolution  of  Ax  =  9.12  km  and  Ay  =  9.38  km  were  performed  to  ensure  that  the 
model’s  behavior  is  robust. 

3.2.4  Nondimensional  Parameters 

The  model  has  the  following  nondimensional  parameters: 

•  the  nondimensional  beta,  defined  as  the  ratio  of  the  advective  time  scale  to  the 
barotropic  Rossby  wave  period:  (3  =  j3dimL2/U ; 
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•  the  inverse  deformation  radius,  defined  as  the  ratio  of  a  horizontal  length  scale 

to  the  internal  baroclinic  deformation  radius:  7  =  y/F[  +  F2  =  L/Ldef ,  where 
the  i— th  layer  Froude  number  is  given  by  F{  =  -  \  ; 

9 

•  the  depth  ratio:  6  =  Di/D2] 

•  the  ratio  of  the  maximum  of  the  lower  layer  inflow  streamfunction,  ^2,  to  the 
maximum  of  the  upper  layer  inflow  streamfunction,  :  a  =  max(^2)/max(^i), 

•  the  nondimensional  interfacial  or  bottom  friction  coefficient,  which  is  the  ratio 
of  the  advective  time  scale  to  the  frictional  decay  time  scale:  A2  =  A2imL/U , 

•  the  nondimensional  biharmonic  diffusion  coefficient,  which  is  the  ratio  of  the 
advective  time  scale  to  the  biharmonic  diffusion  time  scale:  A6  =  AfmL^zU~l. 

3.2.5  Forcing  and  Boundary  Conditions 

The  model  is  forced  by  the  set  of  prescribed  inflow  and  outflow  boundary  condi¬ 
tions.  Such  a  forcing  isolates  baroclinic  and  barotropic  instabilities  as  mechanisms 
responsible  for  low-frequency  variability  and  multiple  dynamical  regimes. 

Two  thin  jets  of  equal  strength  enter  the  domain  from  the  south  and  the  north 
(Figure  3-1),  flow  along  the  closed  western  boundary,  and  collide  in  the  middle  of  the 
domain,  thereby  creating  an  eastward  jet.  The  latter  propagates  across  the  domain 
and  leaves  through  the  opened  eastern  boundary.  The  other  two  regions  of  inflow 
are  located  in  the  northern  and  southern  parts  of  the  eastern  boundary,  where  the 
entering  flow  has  weak  westward  zonal  velocity.  The  geometry  of  the  model  allows  for 
the  interaction  between  the  western  boundary  current  and  the  recirculations  generated 
by  the  instabilities  of  the  eastward  jet. 
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Figure  3-1:  The  geometry  of  the  colliding  jets  model. 


The  inflow  and  the  outflow  on  the  open  boundaries  are  implemented  using  a 
given  streamfunction  field,  which  is  kept  fixed  in  time: 

0  for  {(x,  y)  |  x  =  xw,  ys  <  y  <  yn}  , 

tyeast  for  {(x,  y)  I  X  =  Xe,  Vs  <  y  <  yn i, 

M*,v)  =  Ri{  1  U  J  (3.4) 

*!outh  for  {(xiV)  I  XW  <  X  <  xe,  y  =  ys}  , 

xfysoiith  for  {(^y)  |  xw  <  x  <  xe,  y  =  yn}  , 


V(fcVi(x,y)  =  Rt 


for  {(x,  y)  |  x  =  xw,  ys  <  y  <  yn}  , 
for  {(x,  y)  |  x  =  xe,  ys  <  y  <  yn}  , 


for  {(x,  y)  I  X  =  xe,  ys  <  y  <  yn}  , 
£h*rth  for  {(x,  y)  |  xw  <  x  <  xe,  y  =  ys}  , 
for  {(®>  V)  I  xw  <x  <xe,  y  =  yn }  , 


where  k  =  2, 4,  =  1  and  i?2  —  cr.  The  coefficient  a:  varies  between  zero  and 

one,  thus  allowing  the  creation  of  inflow/outflow  conditions  with  different  degrees  of 
baroclinicity. 
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We  chose  a  free-slip  boundary  condition  for  the  closed  western  boundary. 
Berio ff  and  McWilliams  [3]  compared  the  free-  and  no-slip  boundary  conditions  in  a 
framework  of  a  three-layer  wind-driven  quasi-geostrophic  model  with  lateral  diffusion 
and  showed  that  the  main  difference  in  the  resulting  dynamics  was  that  for  computa¬ 
tionally  achievable  frictional  parameters,  the  western  boundary  current  did  not  have 
local  instability  with  a  free-slip  boundary  condition. 

Here,  the  boundary  conditions  are  chosen  in  such  a  way  that  the  geostrophic 
transport  into  the  domain  balances  the  geostrophic  transport  out  of  the  domain. 
However,  the  total  mass,  which  also  includes  ageostrophic  fields,  is  not  fixed  in  time. 
The  presence  of  open  boundaries  allows  for  small  fluctuations  of  mass  due  to  variations 
in  the  order  of  Rossby  number  fields.  The  description  of  such  fluctuations  is  presented 
in  Chapter  4. 

There  is  zero  potential  vorticity  flux  through  the  boundaries: 


where  iq  —  ( Ui,Vi)T  —  (— ^r)r  is  a  quasi-geostrophic  velocity  vector  and  n  is  a 
vector  normal  to  the  boundary  Cl.  The  contribution  from  planetary  vorticity,  (3y ,  is 
zero  because  the  open  boundaries  are  symmetric:  yn  =  -ys.  The  potential  vorticity, 
V2V>i  +  Fi(ip3-i  -tpi),  also  integrates  to  zero  since  the  streamfunction  on  the  northern 
boundary  is  antisymmetric  to  the  streamfunction  on  the  southern  boundary. 

However,  the  potential  vorticity  flux  into  the  southern  (northern)  part  of  the 
domain  depends  on  the  particular  form  of  the  inflow  and  outflow  conditions.  Different 
conditions  may  result  in  an  intensified  or  a  weakened  potential  vorticity  exchange 
across  the  potential  vorticity  front  associated  with  the  eastward  jet. 

3.2.6  Sponge  Operator 

Traditionally,  models  with  open  boundaries  use  so-called  radiation  condition  schemes 
on  the  outflow  boundaries.  Such  radiation  conditions  attempt  to  allow  the  outgoing 
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part  of  the  disturbance  to  propagate  out  of  the  domain  without  reflecting  back  and 
therefore  without  contaminating  the  interior  solution. 

However,  radiation  conditions  were  not  used  in  the  present  model  for  the  fol¬ 
lowing  reasons.  The  open  eastern  boundary  in  the  model  has  regions  of  both  inflow 
and  outflow  (Figure  3-1).  A  radiation  condition  can  be  applied  only  on  the  outflow 
part  of  the  boundary;  the  inflow  is  fixed  in  time,  thus  providing  steady  forcing.  Nat¬ 
urally,  the  information  on  the  outgoing  disturbance  that  comes  from  the  interior  of 
the  domain  and  leaves  through  the  outflow  part  of  the  boundary  is  independent  from 
the  information  that  enters  the  domain  from  “outside”  through  the  inflow  part  of  the 
boundary.  This  results  in  a  discontinuity  between  the  inflow  and  outflow  values  of 
the  field  to  which  a  radiation  condition  is  applied  (here,  the  streamfunction  field). 

As  a  result,  in  our  model,  the  problem  of  outgoing  radiation  is  resolved  by 
combining  prescribed  inflow/outflow  boundary  conditions  with  an  absorbing  sponge 
layer,  which  is  represented  in  Equation  (3.1)  as 


c  _  a/  dSfdipj  dVj  dS . dipj  dVj 

Si  —  S(x,y)(qi  Qi)  „  (  „  ~  )  o  (a..  a,.)’ 


(3.7) 


dx  dx  dx  dy  dy  dy 
where  S(x,  y)  is  a  sponge  function,  qi  and  ipi  are  instantaneous  potential  vorticity  and 
streamfunction  fields,  and  Qi  and  are  prescribed  “relaxation”  potential  vorticity 
and  streamfunction  fields.  The  latter  are  chosen  as  initial  fields.  Therefore,  in  addi¬ 
tion  to  absorbing  the  outgoing  radiation,  the  sponge  layer  also  allows  for  a  smooth 
transition  of  the  instantaneous  interior  flow  to  the  prescribed  distribution  of  inflow 
and  outflow  boundary  conditions. 

A  particular  form  of  the  sponge  function  is  given  by 


S(x,y) 


U  T, 


decay 


0 

(  x-a  \2  ,  (_U=b_)2 
'Wmt'  ^  y  WnoTth  > 

(  x-a  \2  i  /  V+b  \2 
W™  t>  “r  ^  Wsouth  ) 


for  x  <  a,  —b<y<  6, 
for  x  >  a,  y  >  b, 
for  x  >  a,  y  <  —b, 


(3.8) 


where  Tdecay  =  24  hours  is  the  minimum  decay  scale  in  the  sponge  layer,  a  =  xe  —  30Ax 
and  b  =  yn  —  25A y.  The  corresponding  dimensional  width  of  the  sponge  layer  in  the 
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eastern  part  of  the  domain  is  Weast  =  550  km  and  the  width  in  the  northern  and 
southern  parts  of  the  domain  is  Wnorth  =  Wsouth  —  470  km  (Figure  3-2). 


3.2.7  Numerical  scheme 


The  governing  equations  are  discretized  in  time  with  a  centered  leapfrog  scheme  to 
integrate  forward  the  potential  vorticity  field: 


4?  = 


+  2A£(j7^  +  Si/  -  Vl/), 


where  k  and  l  denote  grid  points  in  x  and  y  directions,  t  —  1,  t ,  and  t  +  1  are 
consecutive  time  steps,  Jh  is  a  Jacobian  term,  Sj.'J1  is  a  sponge  function  term,  and 
V^}1  is  a  dissipation  term.  To  suppress  a  computational  mode  associated  with  the 
leapfrog  scheme,  the  potential  vorticity  fields  at  t  -I- 1  and  t  —  1  are  averaged  every  20 
iterations  to  get  a  corrected  potential  vorticity  field  at  t  +  1. 

To  calculate  Jh,  the  Arakawa  Jacobian  formulation  is  used.  This  allows  no 
net  advection  out  of  the  interior  region  of  potential  vorticity,  enstrophy  or  energy  as 
long  as  the  streamfunction  is  constant  along  the  boundary. 


3.2.8  The  Inversion  of  the  Potential  Vorticity  Field 

The  time  stepping  scheme  is  applied  to  the  potential  vorticity  field  which  is  then 
inverted  to  get  the  streamfunction.  The  inversion  problems  are  solved  separately  for 
the  barotropic,  At  =  F\A  +  F2^i,  and  baroclinic,  V’ftc  =  4>2  ~  Vh)  modes  and  are 
given  by 

V2V>W  =  Fi<72  +  F2qi  -  2y  in  Q,  (3.10) 

i>bt  =  Flfo  +  F2-tpi  in 


and 


(V2  -  7 2)^bc  =  92-9i  in  f2, 
Ac  =  A  in  ft. 


(3.11) 
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The  upper  and  lower  layer  streamfunctions  are  then  restored: 


A  =  (At  +  (~l)lFiAc)l  2-  (3.12) 

FORTRAN  77  routines,  H2GCIS  and  H2GCSS,  from  the  CRAYFISH  package 
were  used  for  solving  the  Helmholtz  Equations  (3.10)  and  (3.11). 


3.2.9  Initialization  of  the  Model 


To  ensure  the  stability  of  currents  entering  the  domain  through  the  open  boundaries, 
the  model  is  initialized  with  a  streamfunction  field  constructed  from  the  Fofonoff 
solution.  The  latter  satisfies  Arnold’s  sufficient  condition  for  stability  of  a  stationary 
nonparallel  flow  [2],  since  >  0. 

To  create  an  initial  condition,  we  first  define  the  following  streamfunction  field 
(Figure  3-3a): 


-1>f  for  x*w  <  x  <  x*e,  0  <  y  <  yn 
>  ipF  for  x*w<x<x*e,  ys  <  y  <  0, 


(3.13) 


where  x*  =  —  x*w  =  L*x/2Ly ,  L*x  =  5000  km  and  tpF  is  the  Fofonoff  solution  satisfying 

=  Py°  [  (y  +  y°)  _  sinh  k(y  +  y0)  +^y  (~l)w  coshfcna;  n?r(y  +  y0)  \ 

k2  \  yo  sinh%0  ^  n  \kn)  coshfcnx0  y0  j 

where  y0  =  x0  -  L*x/2Ly,  k2n  =  k2  +  (n7r/y0)2,  k2  =  y2  +  C,  C  >  0. 

The  potential  vorticity  field,  <7*  =  /3y  +  V2^*  —  y2^*,  has  a  discontinuity  in  the 
middle  of  the  domain  (Figure  3-4),  since 


lim  q(x,y)  =  -  Urn  q(x,y)  ^  0. 

y-+0~  y->0+ 

The  discontinuity  may  produce  instability,  and  therefore  we  smooth  it  out  by  interpo¬ 
lating  q *  between  its  values  at  y  =  -0.1  and  y  =  0.1  to  get  a  new  potential  vorticity 
field,  </**  (Figure  3-3d).  The  inversion  of  this  field  results  in  a  streamfunction, 
(Figure  3-3c)  that  corresponds  to  a  600-km  wide  jet.  The  initial  condition  thus  is 
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Figure  3-4:  Solid  line:  the  meridional  profile  of  potential  vorticity,  qt  =  j3y  + 
72^,,  where  ip*  is  given  in  Equation  3.13.  Dashed  line:  the  meridional  pro 
interpolated  potential  vorticity,  q **. 
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Figure  3-5:  in  the  middle  of  the  domain. 
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chosen  as  ip  =  ip*t  for  (x,  y)  G  £2  (Figure  3-3  e,f).  Notice  that  after  the  interpolation, 
the  necessary  condition  for  instability  is  satisfied,  since  there  are  regions  where 
is  negative  (Figure  3-5).  Smoothed  field  is  no  longer  a  steady  solution. 


3.3  Energetics  of  the  Model 

3.3.1  Energy  Equation 

Multiplying  the  upper  and  the  lower  layer  potential  vorticity  equations  (3.1)  by  —  ipi  jj- 
and  correspondingly  and  adding  them  together  give  the  following  energy  equa- 


+  V-  P  =G. 


(3.14) 


Here,  E  is  the  total  energy  defined  by  the  sum  of  the  depth-integrated  upper  and 
lower  layer  kinetic  energy  and  the  available  potential  energy: 


E  =  tdi  Ut£)  + 


(3.15) 


where  F0  =  FLdi+Fidi ,  ^  %  =  1,2.  The  x—  and  y—  components  of  the  energy 

flux  vector,  P  =  P(x)i  +  P(y)j,  are  defined  by 


P( x)  ^  !  ditpj 


.0, 

Ui  +  g  V’i 


(3.16) 


pm  = 

where  Hi  =  qi  —  (3y.  The  dissipation  term  is  given  by 


(3.17) 


G  —  —  Y.  ditpiiVi  +  Si), 


(3.18) 


where  Vi  and  Si  are  defined  in  section  3.2.1.  In  the  absence  of  dissipation,  the  time 


changes  of  the  total  energy  are  only  due  to  the  divergence  of  the  energy  flux  vector. 
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3.3.2  Energy  Sources  and  Sinks 


Total  energy  sources  and  sinks  can  be  identified  by  integrating  equation  (3.14)  over 
the  model  domain: 

<9E  fxe  ryn 

Q-j.  ~  Jx  2/n)  Vs))  dx  —  J  P(x){xe,y)  dy + 

fXe  ryn 

/  /  G(x,y)  dx  dy,  (3.19) 

J  Xxu  ''Vs 

where  E  =  J^n  E  dx  dy. 

The  first  two  integrals  on  the  right-hand  side  of  equation  (3.19)  are  energy 
sources  provided  by  inflows  through  the  open  boundaries,  whereas  the  last  inte¬ 
gral  is  the  energy  sink  due  to  dissipation.  It  is  important  to  recognize  that  the 
time  dependence  of  the  energy  source  integrals  stems  from  the  following  terms: 
^i{xe,y)%-tVi{xe,y),  ipi(x,yn)-§-tUi(x,yn)  and  ^i{x,ys)f-tUi{x,ys).  Due  to  the  presence 
of  a  strong  sponge  layer,  which  acts  to  damp  perturbations  from  the  prescribed  field, 
'ffi,  the  acceleration  terms  are  negligible.  Therefore,  the  energy  input  due  to  in¬ 
flow/outflow  boundary  conditions  is  time-independent,  which  implies  that  time  vari¬ 
ations  of  total  energy  are  only  due  to  correlation  between  the  solution  and  the  dissi¬ 
pation  term. 

In  contrast,  in  the  wind-driven  quasi-geostrophic  models,  even  though  the 
wind-stress  curl  is  steady,  the  energy  input,  defined  as  the  correlation  between  the 
streamfunction  and  the  wind-stress  curl,  —  fx  fy  ip  •  V  x  r  dx  dy,  can  vary  in  time, 
depending  on  the  solution  structure.  Primeau  [29]  showed  a  remarkable  similarity  be¬ 
tween  the  steady  and  the  time-dependent  solutions  of  a  wind-driven  reduced-gravity 
model  of  McCalpin  and  Haidvogel  [25]  and  estimated  the  differences  in  the  energy 
input  for  multiple  equilibria  at  13%.  Scott  and  Straub  [34]  found  that  antisymmet¬ 
ric  steady  solutions  with  strong  jet  penetration  and  the  jet  axis  aligned  with  the 
maximum  wind  stress  were  the  most  efficient  at  maximizing  the  energy  input.  Non- 
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symmetric  steady  solutions,  on  the  other  hand,  were  less  correlated  with  the  wind 
forcing  and  thus  provided  smaller  energy  input. 

In  their  time-dependent  model,  McCalpin  and  Haidvogel  [25]  found  that  the 
symmetry  of  the  global  wind  forcing  structure  was  important  for  maintaining  the  high 
energy  dynamical  regimes,  which  closely  resembled  the  antisymmetric  steady  state 
solutions  of  Primeau  [29].  For  high  values  of  the  wind  forcing  asymmetry  parameter, 
the  high  energy  state  disappeared,  and  the  flow  was  characterized  by  the  low  energy 
state  with  weakly  penetrating  jet  and  strong  meanders. 

Consider  the  following  thought  experiment  illustrating  the  importance  of  global 
wind  forcing  for  the  existence  and  maintenance  of  multiple  dynamical  regimes.  As¬ 
suming  all  parameters  and  the  basic  state  are  the  same  in  the  wind-driven  and 
boundary-forced  models,  the  linear  vicinities  of  a  fixed  point  corresponding  to  the 
antisymmetric  solution  should  be  topologically  equivalent  for  the  two  models,  since 
the  wind-forcing  drops  out  of  the  linear  equations  of  motion.  As  we  show  later,  for 
the  antisymmetric  steady  solutions,  the  structure  and  the  growth  rate  of  the  most 
unstable  modes  are  indeed  similar  for  the  regional  and  wind-driven  reduced-gravity 
models.  Now,  suppose  the  wind-driven  model  trajectory  is  recovering  from  a  low 
energy  state  and  is  moving  towards  a  high  energy  state,  thus  approaching  a  linear 
neighborhood  of  an  antisymmetric  fixed  point.  In  this  case,  the  recirculations  be¬ 
come  larger,  and  the  jet  straightens.  The  spatial  structure  of  the  solution  is  now 
strongly  correlated  with  the  antisymmetric  wind  stress  curl,  and  thus  the  energy  in¬ 
put  is  increased.  Therefore,  if  the  model  trajectory  is  close  to  the  fixed  point,  but  the 
nonlinear  interactions  are  still  at  work,  the  wind  forcing  could  play  a  crucial  role  in 
further  “directing”  the  model  towards  a  steady  antisymmetric  solution.  The  dissipa¬ 
tion  is  also  important,  and  perhaps  it  might  tend  to  move  the  trajectory  away  from 
the  fixed  point.  Numerical  experiments  showed  that  the  model  initialized  with  the 
antisymmetric  solution  would  take  more  than  20  years  to  move  into  a  lower  energy 
state.  Therefore,  strong  correlation  between  global  wind  forcing  and  an  antisym- 
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metric  solution,  which  expresses  itself  in  the  increased  energy  input,  provides  higher 
chances  of  reaching  and  maintaining  a  high  energy  state  in  a  time-dependent  wind- 
driven  model.  Such  a  mechanism  for  restoring  the  symmetry  of  a  solution  is  absent 
in  the  boundary-forced  model,  where  the  energy  input  in  fixed  in  time.  It  is  only 
nonlinear  eddy-mean  flow  and  eddy-eddy  interactions  and  dissipation  that  affect  the 
model  trajectory  as  it  moves  toward  the  fixed  antisymmetric  point.  We  therefore 
expect  that  in  the  time-dependent  colliding  jets  model,  it  would  be  more  difficult  to 
reach  and,  subsequently,  to  maintain  a  high  energy  state. 

3.3.3  Energy  Convergence  Between  the  Time  Mean  Flow  and 
Perturbations 

The  energy  conversion  rates  between  the  time  mean  flow  and  perturbations  are  impor¬ 
tant  for  understanding  the  dynamics  of  multiple  states  and  of  the  transition  between 
them.  To  derive  the  energy  exchange  rates,  let  us  divide  the  potential  vorticity  and 
streamfunction  fields  into  the  time  mean  (denoted  by  bar)  and  the  perturbation  (de¬ 
noted  by  prime): 

Qi(xty,t)  =  qi(x,y)+q'l(x,y,t),  (3.20) 

ipi(x,y,t)  =  ^i(x,y)+jp'l(x,y,t),  (3.21) 

where  *  =  1,2  are  layer  indices.  Then  the  potential  vorticity  equation  (3.1)  can  be 
rewritten  as 

dq'  _  _  _  _ 

Qf.  Jfyii  Qi)  T  Qi)  +  Jfyii  Qi)  +  JbPii  Qi)  =  'D i  +  T>\  +  S{  +  S'.  (3.22) 

Assuming  that  perturbations  are  small  compared  to  the  time  mean  fields,  the  lin¬ 
earized  perturbation  potential  vorticity  equation  has  the  following  form: 


D+S  D+S' 

Figure  3-6:  General  energy  conversion  diagram  for  the  2-layer  model. 


Multiplying  equation  (3.23)  by  -ip^di  and  summing  over  two  layers  give  the  following 
equation  for  perturbation  energy  density: 


I[|*H  = 


,  cVjx  r/T  \1 

<■  +  ■W*  + 


If  H  §>+^-f>- 
If  Kf)+^f»]+ 

Wi  - V4)  [J'Mt'M  -  J(AW]  - 
Y.d^'M  +  s'd 


(3.24) 


where  K[  =  j  +  J  is  eddy  kinetic  energy  density  and  P'  =  y'  (t)  - C’))- 

is  eddy  potential  energy  density. 

Since  the  perturbation  streamfunction  is  zero  on  the  boundaries,  all  the  diver¬ 
gence  terms  of  equation  (3.24)  disappear  in  the  domain  integral.  Terms  responsible 


101 


for  time  changes  in  the  perturbation  energy  are  the  barotropic  energy  conversion 
rate,  BTC ,  between  the  mean  flow  and  perturbations  due  to  work  by  eddy  momen¬ 
tum  fluxes  and  the  baroclinic  conversion  rate,  BCC,  between  the  mean  flow  and 
perturbation  due  to  work  by  eddy  heat  fluxes: 


BTC(x,y)  = 


V'  j  j/  ,  t  diPi\  .  t  (—  dip'i 

ldy  dy)  +  J{^  dy]  ’ 


(3.25) 


BCC(x,  y)  =  ipo(^-V’2)[^,V'2)-^2>^i)]-  (3.26) 

When  the  conversion  rates  integrated  over  the  model  domain  are  positive,  the  energy 
associated  with  perturbations  grows  in  time,  and  therefore  eddies  draw  the  energy 
from  the  mean  flow.  As  we  noted  earlier,  the  rate  of  external  forcing  due  to  inflow 
boundary  conditions  is  given  by  the  energy  flux  vector  V  •  P,  where  P  is  defined  in 
(3.16)-(3.17).  The  dissipation  rate  for  the  mean  flow  and  for  the  perturbation  are 

2 

-Y,diA{Vi+Si)  (3.27) 

i=  1 

and 

-£  «(»,'  +  S').  (3.28) 

i=l 

Integrals  of  the  above  conversion  rates  over  the  model  domain  are  summarized  in  the 
energy  conversion  diagram  in  Figure  3-6. 


3.4  Ring  Dynamics  at  the  Western  Wall 

The  behavior  of  rings  and  their  interaction  with  the  western  boundary  currents  may 
be  crucial  the  stability  of  the  system  and  therefore  on  the  transition  between  multiple 
dynamical  states. 
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Figure  3-7:  A  schematic  illustrating  the  behavior  of  rings  on  the  western  boundary. 

A  ring  trajectory  at  the  western  boundary  is  affected  by  two  competing  pro¬ 
cesses  (Figure  3-7).  The  first  one  is  the  formation  of  an  image  vortex  with  the  opposite 
vorticity  sign  due  to  a  rigid  wall  boundary  condition.  A  newly  created  dipole  struc¬ 
ture  tends  to  propel  itself  along  the  boundary  northward  for  anticyclonic  rings  and 
southward  for  cyclonic  rings.  The  second  process  affecting  the  ring  movement  is  its 
advection  towards  the  middle  of  the  domain  by  the  western  boundary  current.  There¬ 
fore,  depending  of  the  strength  of  the  ring,  it  either  moves  towards  the  eastward  jet 
or  away  from  it. 

If  the  ring  is  weak  enough  and  therefore  is  advected  towards  the  middle  of  the 
domain,  it  usually  rejoins  the  eastward  jet  and  may,  by  injecting  anomalous  potential 
vorticity,  further  perturb  it.  However,  if  the  ring  is  strong  enough  to  overcome  the 
advective  force  of  the  western  boundary  current,  two  scenarios  are  possible.  First, 
the  ring-image  vortex  dipole  moves  along  the  wall  until  entering  an  absorbing  sponge 
layer.  Such  a  ring  might  be  able  either  to  halt  or  to  change  the  path  of  the  western 
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boundary  current.  This,  in  turn,  might  break  the  symmetry  of  the  inflow  conditions 
near  the  separation  point  and  lead  to  the  overshoot  of  the  opposite  flowing  boundary 
current.  Second,  a  weaker  ring  may  dissipate  enough  energy  while  traveling  along  the 
wall,  to  reverse  its  trajectory  before  reaching  the  sponge  layer  and  subsequently  to 
join  the  eastward  jet. 
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Chapter  4 


Results  for  Reduced- Gravity 
Colliding  Jets  Model 

4.1  Model  Formulation 

4.1.1  Governing  Equation 

This  chapter  reports  the  results  for  the  reduced-gravity  version  of  the  colliding  jet 
model  which  isolates  the  importance  of  barotropic  instability  and  dissipation  for 
the  existence  and  maintenance  of  low  frequency  variability  and  multiple  dynamical 
regimes  of  a  quasi-geostrophic  jet  with  recirculations.  The  governing  potential  vor- 
ticity  equation  is  easily  deduced  from  equation  (3.1)  when  the  lower  layer  is  assumed 
infinitely  deep  and  motionless: 

||  +J(M=V  +  S  in!).  (4.1) 

subject  to  the  boundary  conditions: 


Here,  ip  —  ip(x,y )  is  a  quasi-geostrophic  streamfunction,  q  =  V2ip  +  /3y  —  7 2ip  is 
potential  vorticity.  The  parameter  space  is  reduced  to  , 


L2  „  pdimL 2  „ 


AfmL  Afrn 
U  ’  ^  ~ 


since  the  depth  ratio,  <5,  and  the  baroclinic  coefficient,  a,  are  zero.  Coefficient  .4.2 
represents  the  interfacial  friction,  and  therefore  the  dissipation  in  the  upper  layer 
is  given  by  V  —  -A^ip  -  A&V^ip.  The  model  geometry  and  the  structure  of  the 
sponge  layer  are  as  described  in  Chapter  3. 

4.1.2  Energetics 

For  the  case  of  infinitely  deep  and  motionless  lower  layer,  the  energy  equation  is 

dE 

aF  +  v '  p  =G-  («) 

Here,  E  is  the  total  energy  defined  by  the  sum  of  the  kinetic  and  available  potential 


energy: 


1  (f  di‘  \  2  / di’Y  2  ,2 

E  =  2  H&J  W  +7V’ 


The  x—  and  y—  components  of  the  energy  flux  vector,  P  =  P^i  +  P(y)j,  are  defined 


p{x)  =  -  i>Tlu  —  ~ip2 


W  Uj 

P(y)  =  ^~Qt  ”  ^ 

where  n  =  q  —  (3y.  The  dissipation  term  is  given  by 


G  =  -if>{V  +  S). 
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Figure  4-1:  Energy  diagram  for  the  reduced-gravity  model. 


Since  there  is  no  horizontal  heat  transfer,  the  baroclinic  conversion  term  between  the 
perturbation  and  the  time  mean  state  is  zero.  It  is  only  due  to  work  by  Reynolds 
stresses  that  there  is  an  energy  exchange  between  the  mean  flow  and  perturbations, 
given  by  the  barotropic  conversion  rate: 


The  energy  input  due  to  inflow  boundary  conditions  is  fixed  in  time  and  is  given  by 
domain-integrated  divergence  of  the  energy  flux  vector,  P.  The  dissipation  rate  for  the 
mean  flow  and  for  perturbations  are  —^{T>  +  S )  and  —ip'(V  +  S')  respectively.  The 
energy  conversion  diagram  (Figure  4-1)  summarizes  the  domain-integrated  energy 
transfer  rates  for  the  reduced-gravity  model. 


4.2  Linear  Stability  Analysis  of  Steady  Antisym¬ 
metric  Solutions 


4.2.1  Steady  States  vs  Interfacial  Friction  Coefficient 


In  this  section,  we  describe  steady  antisymmetric  solutions  of  the  reduced-gravity 
model.  As  we  show  later,  such  solutions  closely  resemble  variability  associated  with 
a  high  energy  state  of  the  time-dependent  model.  Therefore,  the  linear  stability 
characteristics  of  these  solutions  are  important  for  understanding  the  eddy-mean  flow 
interaction  present  in  a  high  energy  state. 

Steady  antisymmetric  solutions  of  the  reduced-gravity  model  were  calculated 
by  adapting  a  version  of  the  numerical  code  with  enforced  meridional  antisymmetry. 
On  each  iteration,  the  potential  vorticity  field  was  stepped  forward  in  time  only  for 
the  southern  half  of  the  domain,  y  <  0,  subject  to  if>  =  0  at  y  =  0.  The  fields  for 
the  northern  half  of  the  domain,  y  >  0,  were  taken  as  mirror  images  of  the  fields 
f°r  y  <  0.  The  model  was  integrated  until  the  following  numerical  equivalent  of  an 
equilibration  criteria  was  satisfied: 


\E(t  +  At)  -  E{t)\ 

m 


<  1(T9A  t, 


where  E(t )  is  total  energy  at  time  t  and  At  is  the  time  step. 

Figure  4-2  presents  steady  antisymmetric  solutions  of  the  reduced-gravity  model 
for  Afm  =  8 -1010  m4  sec-1,  Tin  =  16  Sv  and  the  following  values  of  the  interfacial  fric¬ 
tion  coefficient,  Afm:  l.O-KT7,  5.0-10"8,  3.5-KT8, 3.0-lCr8, 2.8-HT8, 2.5-KT8  sec"1. 


With  an  increase  in  the  dimensional  interfacial  friction  coefficient,  Afm,  the  total  en¬ 
ergy  level  of  the  steady  solutions  diminishes  from  0.08  to  0.052  (Figure  4-3a),  which 
is  directly  related  to  a  decrease  from  2350  to  1000  km  in  the  longitudinal  size  of  the 
recirculating  gyres  (Figure  4-3b).  The  maximum  in  the  total  eastward  (westward) 
transport  also  drops  from  41  to  26  (31  to  16)  Sv  (Figure  4-3c,d),  indicating  weaker 
eastward  jet  and  recirculations. 
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Figure  4-2:  Streamfunction  contours  of  steady  antisymmetric  solutions  for  Afm  = 
1.0  •  10-7  sec-1  (a),  5.0  •  10~8  sec-1  (b),  3.5  •  10~8  sec-1  (c),  3.0  •  10-8  sec-1  (d) 
2.8  •  10~8  sec-1  (e),  2.5  •  10-8  sec-1  (f).  Parameters:  T{n  =  16  Sv,  Afm  =  8  •  1010 
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Figure  4-3:  Total  nondimensional  energy  (a),  the  longitudinal  size  of  the  recirculations 
(b),  and  maxima  in  the  total  eastward  (c)  and  westward  (d)  transport  of  steady 
antisymmetric  solutions,  ^5,  of  the  reduced-gravity  model  as  functions  of  dimensional 
interfacial  drag  coefficient,  Afm.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1. 


4.2.2  Growth  Rates,  Phase  Speeds  and  Structures 


of  Unstable  Modes 


To  determine  the  growth  rates,  phase  speeds  and  spatial  structures  of  growing  modes 
associated  with  steady  state  solutions,  we  used  an  initial  value  technique  applied  to  the 
original  nonlinear  model.  A  similar  approach  was  adopted  by  Simmons  and  Hoskins 
[36]  in  a  study  of  quasi-geostrophic  baroclinic  instabilities  on  a  sphere.  The  model  was 
initialized  with  the  sum  of  a  steady  state  solution,  fl/s,  and  a  uniformly  distributed 
random  noise,  'ip'ix,  y).  The  amplitude  of  the  perturbation  potential  vorticity  field 
was  chosen  as  1%  of  the  maximum  in  the  steady  state  potential  vorticity  field,  Qs- 
Let  us  represent  the  perturbation  streamfunction  as  a  modal  solution  of  the 
form  ip'(x,y,t)  =  Real  A(x,y)e(~ar+l<T^t,  where  A(x,y)  =  AT(x,y)  +iAi(x,y)  is  a 
spatially  varying  amplitude,  Ar  and  A;  are  real  functions,  oT  is  a  growth  rate  and  eq 
is  frequency.  The  total  energy  of  the  perturbation  field  is  given  by 

E(i,'(x,y,t))  =  ^  +  72|B|2},  (4.9) 

where  B  =  AT  cos  crtf  —  Ai  sin  a  it  +  i  (A,  coscqf  +  Ar  sino*) .  Assuming  that  the 
perturbation  amplitude  is  small,  and  therefore  the  nature  of  the  eddy-mean  flow 
interaction  is  linear,  the  growth  rate  can  be  estimated  by  fitting  a  straight  line  to  the 
natural  logarithm  of  the  perturbation  energy: 


{r)B  2  SB  2  1 

^  +12\B\2) 


(4.10) 


When  Afm  =  1.0- 10-7  sec-1,  the  perturbation  energy  decays  to  zero  with  time, 
and  therefore  the  corresponding  steady  antisymmetric  solution  is  stable.  However, 
for  Afm  >  5.0- 10-8  sec-1,  steady  antisymmetric  solutions  are  unstable,  and  the  time 
evolution  of  \nE(xJ;'(x,  y))  (Figure  4-4a)  is  characterized  by  three  distinct  regions 
which  reflect  changes  in  the  nature  of  internal  interactions.  Initially,  the  uniformly 
distributed  perturbation  field  does  not  correspond  to  a  preferred  mode  of  instability 
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(a) 


(b) 


Figure  4-4:  (a):  Natural  logarithm  of  energy  of  ip'  =  ^  as  a  function  of  time, 
where  Vi's  is  an  antisymmetric  steady  solution  and  is  the  solution  of  the  time- 
dependent  model  initialized  with  a  perturbed  steady  solution.  Interfacial  friction 
coefficient:  A2  =  5.0  •  10-8  sec-1  (thin  solid  line),  3.5  •  10-8  sec'1  (dotted  line), 
3.0  •  10-8  sec-1  (thick  solid  line),  2.8  •  10-8  sec-1  (dashed-dotted  line),  2.5  •  10-8  sec-1 
(dashed  line).  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1,  (b):  Growth  rate, 
<jr,  (in  year-1  )  as  a  function  of  the  dimensional  interfacial  drag  coefficient,  Afm. 
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and  therefore  is  not  efficient  at  extracting  energy  from  the  mean  flow.  After  3  or 
4  years,  the  exponentially  growing  mode  arises  on  the  axis  of  the  eastward  jet  and 
in  the  area  of  recirculations,  and  hence  In  E(tp'(x,  y,  t))  displays  linear  growth.  For 
A% im  <  3.5  •  10-8  sec-1,  the  perturbation  amplitude  becomes  so  large  after  13  to  15 
years,  that  the  nonlinear  eddy-eddy  interactions  become  dominant,  and  the  growth 
rate  saturates.  Notice  that  for  Afm  —  5.0  •  10-8  sec-1,  the  growth  rate  was  still  linear 
after  27  years.  This  indicates  that  a  trajectory  of  the  time-dependent  model  might 
stay  in  the  vicinity  of  the  fixed  point  associated  with  the  antisymmetric  solution  for  a 
long  time.  In  general,  there  might  be  several  unstable,  exponentially  growing  modes 
in  the  system,  which  would  then  enforce  curvature  of  In  E(ip'(x,  y ,  t ))  during  a  stage  of 
exponential  growth.  The  absence  of  such  curvature  in  In  E(ip'(x,y,t))  (Figure  4-4a) 
in  the  present  model  indicates  the  existence  of  just  one  unstable  mode. 

The  e-folding  time,  1/cr,  of  the  unstable,  oscillating  mode  decreases  from  4.8 
years  for  A =  5.0  •  10-8  sec-1  to  417  days  for  A%tm  =  2.5  •  10-8  sec-1  (Fig¬ 
ure  4-4b).  As  we  show  later,  the  time-dependent  regional  model  displays  multiple 
dynamical  regimes  when  Afm  =  3  •  10-8.  The  corresponding  e-folding  time  for  the 
nearly  antisymmetric  solution  is  481  days.  In  comparison,  Primeau  [29]  estimated  the 
e-folding  time  at  897  days  for  the  antisymmetric  steady  state  solution  of  the  wind- 
driven  reduced-gravity  model  with  the  parameter  set  from  McCalpin  and  Haidvogel’s 
model  [25].  The  differences  in  the  e-folding  times  between  the  present  model  and  the 
wind-driven  case  could  be  attributed  to  the  differences  in  the  strength  of  forcing. 

Snapshots  of  the  most  unstable  modes  (Figure  4-5)  show  a  spatial  structure 
similar  to  that  calculated  by  Primeau  [29].  There  are  four  cells  on  the  axis  of  the 
eastward  jet  and  weaker  cells  in  the  area  of  westward  recirculations.  Visualization 
of  the  mode’s  time  series  reveals  a  very  complicated  propagation  pattern,  when  am¬ 
plitudes  and  spatial  extents  of  the  cells  located  on  the  jet  axis  slowly  grow  at  some 
times  and  quickly  shift  westward  at  other  times.  The  westward  propagation  speed  of 
the  cells  in  the  area  of  the  recirculations  is  constant  in  time.  Qualitative  similarity 


113 


Y(  KM ) 


between  the  results  of  linear  stability  analysis  of  the  antisymmetric  solutions  for  the 
regional  and  the  wind-driven  model  stems  from  the  fact  that  wind-forcing  drops  out 
of  the  linearized  equations  of  motion  for  the  wind-driven  case. 

When  nonlinearities  become  large,  the  potential  vorticity  anomalies  associated 
with  four  cells  start  interacting  with  each  other  leading  to  saturation  in  the  energy 
growth  and  to  formation  of  rings  in  the  area  east  of  the  recirculations.  This  process 
might  be  identified  as  a  spatial  growth  of  instabilities. 

4.3  Sensitivity  Study  for  Various  Interfacial  Fric¬ 
tion  Coefficients 

4.3.1  Classification  of  States 

The  regional  colliding  jet  model  showed  strong  dependence  of  the  flow  character  on 
the  value  of  the  interfacial  friction  coefficient,  Afm.  Figures  4-6  and  4-7  summarize 
the  results  of  the  sensitivity  study  for  the  following  two  values  of  the  biharmonic 
diffusion  coefficient:  Afm  =  8  •  1010  m4  sec-1  and  Afm  =  8  •  108  m4  sec-1.  For  both 
values  of  Afm,  with  a  decrease  in  A%im,  the  model  progressed  through  the  following 
regimes: 

•  State  A.  The  flow  structure  consists  of  a  very  broad,  stable  eastward  jet;  the 
recirculations  are  absent. 

•  State  B.  The  flow  structure  is  stable  and  consists  of  an  eastward  jet  surrounded 
by  two  recirculating  gyres. 

•  State  C.  The  flow  structure  is  unstable  and  consists  of  an  eastward  jet  sur¬ 
rounded  by  two  recirculating  gyres.  The  eddy  field  is  dominated  by  eastward 
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A  dim 

STATE 

FLOW  CHARACTER 

5.0  •  10-9 

F 

one  dynamical  state:  eastward  jet  with  large-amplitude, 
westward  propagating  meanders  and  large  meridional 
shifts  in  the  jet  separation  point;  strong  rings  are  formed 
on  the  western  boundary;  recirculations  are  absent 

2.0 -10-8 

F 

as  above 

4.0  •  10-8 

E 

low  frequency  variability  associated  with  two  dynami¬ 
cal  states:  high  energy  state  is  characterized  by  large 
/  medium  /  small  recirculations  and  weak  /  medium  / 
strong  eddy  and  ring  fields,  rings  are  formed  east  of  the 
recirculations;  low  energy  state  has  no  recirculations  and 
is  characterized  by  large-amplitude,  westward  propagat¬ 
ing  meanders,  large  meridional  shifts  in  the  jet  separation 
point  and  strong  rings  formed  on  the  western  boundary 

5.0  •  10-8 

E 

as  above 

6.0  •  10-8 

E 

as  above 

7.5  •  10-8 

D 

low  frequency  variability  associated  with  eastward  jet 
and  longitudinally  pulsating  recirculations,  strong  ring 
formation,  eastward  and  westward  propagating  mean¬ 
ders,  weak  and  medium  rings  are  formed  east  of  the 
recirculations 

1.0- 10-7 

C 

one  dynamical  state:  eastward  jet  and  two  recirculating 
gyres,  weak  ring  formation,  eastward  and  westward  prop¬ 
agating  meanders 

h- 

1 

O 

r-H 

LO 

t“H 

B 

stable  state:  eastward  jet  with  two  recirculating  gyres 

2.0- 10-7 

B 

as  above 

2.5  •  10-7 

A 

stable  state:  broad  stable  jet,  no  recirculations 

Figure  4-6:  Flow  character  vs  interfacial  friction  coefficient,  Afm  (in  sec-1),  when 
Tin  =  16  Sv  and  Afm  =  8  •  108  m4  sec-1. 
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Figure  4-7:  Flow  character  vs  interfacial  friction  coefficient,  Afm  (in  sec  1),  when 
Tin  —  16  Sv  and  Afm  =  8  •  1010  m4  sec-1. 
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and  westward  propagating  meanders.  Weak  rings  form  east  of  the  recirculating 
gyres. 

•  State  D.  The  flow  structure  is  unstable  and  consists  of  an  eastward  jet  and 
two  recirculating  gyres.  The  longitudinal  size  of  the  recirculations  changes  with 
time  thus  giving  rise  to  low  frequency  variability.  The  amplitude  and  the  spatial 
structure  of  the  eddy  field  are  strongly  time-dependent.  The  eddy  field  is  char¬ 
acterized  by  eastward  and  westward  propagating  meanders.  Their  interactions 
give  birth  to  rings  in  the  area  east  of  the  recirculation  gyres. 

•  State  E.  The  flow  structure  displays  well-developed  low  frequency  variability. 
The  model  stays  in  one  of  the  two  dynamical  regimes:  high  energy  state  D  or 
low  energy  state  F. 

•  State  F.  The  flow  structure  is  unstable  and  consists  of  an  eastward  jet  with  large 
meridional  excursions  in  the  separation  point.  The  eddy  field  is  dominated  by 
large-amplitude,  westward  propagating  meanders  and  strong  rings  formed  only 
on  the  western  boundary.  The  recirculations  are  absent  in  instantaneous  fields. 

The  existence  of  multiple  energy  regimes  is  illustrated  in  Figure  4-8,  which 
shows  normalized  histograms  of  total  energy  for  various  values  of  the  interfacial  fric¬ 
tion  coefficient.  The  energy  level  of  steady  antisymmetric  solutions  discussed  in  sec¬ 
tion  4.2.1  is  denoted  by  a  dashed  line  and  is  either  equal  or  higher  than  that  of 
the  corresponding  time-dependent  solutions.  For  weakly  unstable  flows  ( Afm  = 
5  •  10  7  sec  J),  the  trajectory  of  the  time-dependent  model  stays  in  the  vicinity  on 
a  fixed  antisymmetric  point,  and  therefore  the  range  of  diagnosed  energy  values  is 
small  and  is  close  to  E(\ ts)-  The  corresponding  energy  histogram  has  only  one  peak, 
indicating  a  unique  dynamical  regime.  For  stronger  unstable  flows,  the  range  of  en¬ 
ergy  values  of  the  time-dependent  model  becomes  larger.  With  further  decrease  in 
the  energy  histogram  changes  its  shape  from  a  normal  to  a  bimodal  distribu¬ 
tion.  The  model  now  has  two  preferred  energy  states,  each  one  associated  with  a 
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Figure  4-8:  Total  nondimensional  energy  level,  as  diagnosed  from  the  time-dependent 
model,  vs  interfacial  friction  coefficient,  Afm.  Normalized  energy  histograms  are 
shown  on  top  of  each  energy  distribution.  The  dashed  line  denotes  the  energy  level 
of  the  antisymmetric  steady  solutions. 
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distinct  dynamical  regime.  For  small  values  of  the  interfacial  friction  coefficient,  the 
high  energy  state  becomes  unstable,  and  the  model  trajectory  stays  in  a  dynamical 
regime  associated  with  a  low  energy  state. 

4.3.2  Examples  of  States 

In  this  section,  wre  briefly  describe  the  dynamics  of  different  states.  The  full  analysis 
of  the  model  will  be  given  in  section  4.4  for  the  case  of  state  E. 

State  A 

For  high  values  of  the  interfacial  drag  coefficient,  the  flow  has  the  following  structure. 
Part  of  the  meridional  inflow  separates  from  the  western  boundary  current  and  flows 
eastward  creating  a  broad  stable  current  that  leaves  the  domain  though  the  eastern 
boundary.  Another  part  of  the  flow  represents  two  colliding  currents  that  form  a 
narrow  eastward  jet  that  broadens  400-700  km  east  of  the  separation  point  (Figure 
4-9).  There  are  no  recirculations  in  this  state,  although  there  are  weak  (0.2  m  sec-1) 
westward  velocities  north  and  south  of  the  narrow  part  of  the  jet.  Instead  of  closing 
on  themselves,  contours  of  streamfunction  and  potential  vorticity  loop  towards  the 
western  boundary  current  and  then  curve  back  eastward.  Even  though  is  negative 
in  the  eastward  jet,  therefore  satisfying  the  necessary  condition  for  instability,  the 
interfacial  drag  is  so  strong  that  it  suppresses  any  growing  disturbances. 

State  B 

An  example  of  the  flow  structure  for  state  B  is  presented  in  section  4.2.1  (Figure  4-2a) 
for  Afm  =  1  •  10-7sec-1  and  Af™  =  8  •  1010  m4  sec-1. 
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Figure  4-9:  Contours  of  time  mean  streamfunction  (a)  and  potential  vorticity  (b)  for 
A%im  —  2.5  •  10-7  sec-1,  initial  inflow  transport  of  16  Sv  and  Afm  —  8  •  108  m4  sec-1. 
Solid  line  is  used  for  positive  values,  dashed  line  is  used  for  negative  values. 
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Figure  4-10:  Contours  of  time  mean  streamfunction  (a)  and  potential  vorticity  (b) 
for  Ai  —  5.0  •  10-8  sec-1,  initial  inflow  transport  of  16  Sv  and  Ae  =  8  ■  1010  m4  sec-1. 
Solid  line  is  used  for  positive  values,  dashed  line  is  used  for  negative  values. 
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Figure  4-11:  Energy  diagram  for  A2  =  5.0  -10  8  sec  x,  initial  inflow  transport  of  16 
Sv  and  ^4.6  =  8  *  1010  m4  sec-1. 

State  C 

In  agreement  with  linear  stability  analysis  of  section  4.2.1,  the  flow  becomes  unstable 
when  the  interfacial  friction  coefficient  falls  below  some  critical  value.  In  state  C,  when 
Afm  is  close  to  the  critical  value,  the  model  trajectory  is  trapped  in  the  vicinity  of 
an  antisymmetric  fixed  point,  since  the  eddy-mean  flow  interactions  dominate  the 
eddy-eddy  interactions.  This  explains  why  the  model  stays  in  one  dynamical  regime 
and  why  the  time  mean  streamfunction  and  potential  vorticity  fields  (Figure  4-10) 
closely  resemble  the  steady  antisymmetric  solution.  The  flow  structure  consists  of  an 
eastward  flowing  jet  and  two  recirculating  gyres.  The  longitudinal  size  of  the  time 
mean  recirculations  equals  1500  km  (compared  to  1600  km  in  the  case  of  a  steady 
solution)  and  does  not  substantially  change  in  time. 

The  energy  conversion  diagram  summarizing  the  energetics  of  the  model  shows 
that  the  barotropic  conversion  term  is  small.  The  perturbation  energy  is  only  3.3%  of 
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Figure  4-12:  Time  mean  perturbation  kinetic  energy  density  (a),  perturbation  po¬ 
tential  energy  density  (b),  barotropic  conversion  energy  density  (c)  and  perturbation 
energy  dissipation  (d).  Parameters:  Tin  =  16  Sv,  A2  =  5.0  •  10"8  sec-1,  A6  =  8  •  1010 
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the  total  energy.  The  eddy  field  is  dominated  by  eastward  and  westward  propagating 
meanders.  Perturbation  kinetic  and  potential  energy  densities  (Figure  4-12)  have  spa¬ 
tial  structures  resembling  the  structure  of  the  linearly  unstable  mode.  Apparently, 
the  dynamics  is  dominated  by  the  linear  eddy-mean  flow  interactions  and  weakly- 
nonlinear  eddy-eddy  interactions.  The  amplitude  of  the  perturbation  fields  near  the 
separation  point  shows  that  the  latter  has  weak  meridional  excursions  due  to  me¬ 
ander  movement.  The  dissipation  of  the  perturbation  energy  occurs  in  the  area  of 
western  boundary  current  and  on  the  axis  of  the  jet  west  of  1000  km.  The  barotropic 
conversion  energy  density  is  aligned  with  the  jet  axis  and  is  maximum  west  of  1500 
km.  The  spatial  structure  of  the  perturbation  fields  indicates  that  rings  formed  east 
of  the  recirculations  are  so  weak  that  they  are  quickly  dissipated  during  their  west¬ 
ward  propagation  in  the  recirculation  regions.  Therefore,  they  mostly  provide  weak 
potential  vorticity  mixing  across  the  front  in  the  region  east  of  the  recirculations. 

State  D 

In  state  D,  the  model  exhibits  low  frequency  variability  associated  with  changes  in  the 
zonal  extent  of  the  recirculating  gyres.  The  time  mean  streamfunction  and  potential 
vorticity  fields  (Figure  4-13)  have  1200  km  long  recirculations.  The  energy  level  of  the 
time  mean  streamfunction  is  14%  smaller  than  the  energy  level  of  the  antisymmetric 
steady  solution.  During  the  flow  evolution,  the  recirculation  size  fluctuates  between 
700  and  almost  1800  km.  Since  most  of  the  potential  energy  is  contained  in  the 
recirculation  areas,  and  because  the  total  energy  is  dominated  by  the  potential  energy, 
the  time  series  of  total  energy  shows  periods  of  high  and  low  values  (Figure  4-14). 

For  high  values  of  total  energy,  the  model  trajectory  is  in  the  vicinity  of  an  an¬ 
tisymmetric  steady  solution.  The  recirculation  gyres  are  about  1800  km  long,  and  the 
eddy  field  is  very  weak.  The  flow  however,  is  linearly  unstable,  and  therefore  growing 
disturbances  draw  energy  from  the  mean  flow.  The  amplitude  of  the  eddy  field  in¬ 
creases  with  time,  until  nonlinear  interactions  start  dominating.  At  this  moment,  the 
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Figure  4-13:  Contours  of  time  mean  streamfunction  (a)  and  potential  vorticity  (b) 
for  A2  =  3.5  •  10-8  sec-1,  initial  inflow  transport  of  16  Sv  and  A§  =  8  •  1010  m4  sec-1. 
Solid  line  is  used  for  positive  values,  dashed  line  is  used  for  negative  values. 
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Figure  4-14:  Time  series  of  total  nondimensional  energy  for  Ai  =  3.5  -10  8  sec  1, 
initial  inflow  transport  of  16  Sv  and  As  =  8  •  1010  m4  sec-1. 

meanders  reach  large  amplitudes,  which  leads  to  the  formation  of  stronger  rings  east 
of  the  recirculations.  As  the  rings  travel  westward  in  the  recirculation  region,  they 
mix  opposite  sign  vorticity  into  the  surrounding  fluid  and  homogenize  the  surround¬ 
ing,  thus  decreasing  the  recirculation  strength.  Upon  reaching  the  western  wall,  they 
are  advected  by  the  western  boundary  current  and  entrained  into  the  eastward  jet. 
The  formation  and  propagation  of  strong  meanders  and  rings  leads  to  a  decrease  in 
the  size  of  the  recirculating  gyres  and  to  a  subsequent  drop  in  the  total  energy  level. 
Once  the  recirculations  shrink,  the  flow  becomes  less  turbulent,  and  perturbations 
decay.  At  this  point,  the  advection  of  potential  vorticity  by  the  western  boundary 
current  becomes  dominant  and  restores  large  recirculations,  which  causes  the  total 
energy  to  increase. 

The  conical  shape  of  perturbation  kinetic  and  potential  density  fields  (Figure 
4-15)  is  connected  to  the  existence  of  meanders  and  rings  not  only  on  the  axis  of  the 
jet,  but  also  in  the  recirculation  regions  and  along  the  western  wall.  The  barotropic 
conversion  density  and  the  dissipation  of  perturbation  energy  are  maximum  on  the 
axis  of  the  jet  west  of  500  km  and  along  the  western  wall. 
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Figure  4-15:  Time  mean  perturbation  kinetic  energy  density  (a),  perturbation  po¬ 
tential  energy  density  (b),  barotropic  conversion  energy  density  (c)  and  perturbation 
energy  dissipation  (d).  Parameters:  Tin  =  16  Sv,  A2  =  3.5  •  1(T8  sec"1,  A6=8- 1010 
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Figure  4-16:  Energy  diagram  for  A2  =  3.5  -10  8  sec  1,  initial  inflow  transport  of  16 
Sv  and  Ae  =  8  ■  1010  m4  sec-1. 

The  energy  diagram  (Figure  4-16)  summarizes  the  domain-integrated  energy 
transfers  for  state  D.  Compared  to  state  C,  the  energy  associated  with  the  mean 
field  did  not  change  much.  However,  the  perturbation  kinetic  energy  was  an  order  of 
magnitude  larger,  and  the  perturbation  potential  energy  grew  more  than  three  times. 
The  barotropic  conversion  was  almost  twice  as  large.  Since  the  interfacial  friction 
was  weaker,  both  mean  and  perturbation  energy  dissipation  terms  were  smaller  than 
in  state  C.  The  partition  of  perturbation  energy  reached  9%  compared  to  only  3%  in 
state  C. 

State  E 

In  state  E,  the  model  exhibits  strong  low  frequency  variability  associated  with  high 
energy  periods,  when  the  flow  structure  is  as  in  state  D,  and  with  low  energy  periods, 
when  the  flow  structure  is  as  in  state  F.  The  details  of  the  dynamics  for  state  E  are 
presented  in  section  4.4. 
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State  F 


For  small  values  of  the  interfacial  friction  coefficient,  the  flow  is  strongly  unstable. 
After  perturbing  an  antisymmetric  steady  solution  (see  section  4.2.1),  the  eddy  field 
grows  rapidly,  drastically  changing  the  structure  of  the  mean  flow  and  causing  the 
recirculations  to  disappear.  The  separation  point  of  the  eastward  jet  suffers  large 
meridional  excursions.  Since  northward  and  southward  overshoots  are  equally  possi¬ 
ble,  the  western  boundary  currents  appear  to  be  separating  early  from  the  wall  and 
flowing  eastward,  thus  creating  a  region  of  zero  streamfunction  and  “homogenized” 
potential  vorticity  (Figure  4-17).  The  time  mean  fields  are  therefore  misleading  and 
can  not  be  interpreted  as  structures  persisting  in  time.  The  recirculations  are  not 
present  in  the  instantaneous  fields.  Instead,  there  are  individual  vortices,  trapped 
between  meander’s  highs  and  lows,  which  provide  weak  mixing  on  the  flanks  of  the 
jet. 

The  strong  eddy  field  is  dominated  by  large  amplitude,  westward  propagating 
meanders,  which  do  not  allow  any  ring  formation  in  the  interior  of  the  model  domain. 
However,  the  interaction  of  these  meanders  with  the  western  wall  creates  large  excur¬ 
sions  in  the  separation  point  and  strong  rings  that,  due  to  a  no-normal-flow  boundary 
condition,  move  along  the  western  boundary  against  the  advective  force  of  the  west¬ 
ern  boundary  current.  Boundary  ring  formation  provides  strong  potential  vorticity 
mixing  across  the  front  and  tends  to  decrease  the  asymmetry  in  the  inflow  conditions 
near  the  collision  point.  However,  for  a  given  value  of  the  interfacial  friction  coeffi¬ 
cient,  the  meander  amplitudes  are  large  enough  to  guarantee  strong  excursions  in  the 
separation  point,  which  keeps  the  model  trajectory  in  a  low  energy  state. 

In  contrast  with  states  C  and  D,  the  spatial  structure  of  perturbation  energy, 
its  dissipation  and  barotropic  conversion  (Figure  4-18)  have  maxima  along  the  western 
boundary.  Their  penetration  into  the  model  interior  is  less  than  500  km.  The  energy 
diagram  (Figure  4-19)  summarizes  domain-integrated  energy  transfers. 
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Figure  4-17:  Contours  of  time  mean  streamfunction  (a)  and  potential  vorticity  (b) 
for  A2  =  2.0  •  10-8  sec-1,  initial  inflow  transport  of  16  Sv  and  A%  =  8  •  1010  m4  sec-1. 
Solid  line  is  used  for  positive  values,  dashed  line  is  used  for  negative  values. 
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Figure  4-18:  Time  mean  perturbation  kinetic  energy  density  (a),  perturbation  po¬ 
tential  energy  density  (b) ,  barotropic  conversion  energy  density  (c)  and  perturbation 
energy  dissipation  (d).  Parameters:  Tin  =  16  Sv,  A2  =  2.0  •  10-8  sec-1,  A$  =  8  •  1010 

4  —1 

nr  sec  \ 


132 


inflow 


1.4-10 

3 

K  1.3-1 03 

k"  5.1-1 0‘3 

BTC 

1.1-1  o'3 

P  5.T10 

P  1.2-10 

-3.3*1  O'4 

-1.1-ltf 

D+S  D+S' 

Figure  4-19:  Energy  diagram  for  A2  —  2.0  •  10-8  sec-1,  initial  inflow  transport  of  16 
Sv  and  A&  =  8  •  1010  m4  sec-1. 

4.4  Reference  Run:  State  E 

4.4.1  Dimensional  Parameters 

As  a  reference  case,  we  chose  a  run  that  displayed  well-developed  multiple  dynamical 
regimes.  The  dimensional  scales  and  parameters  were  chosen  as  follows: 

L  =  2400  km, 

Lx  —  4667  km, 

Ldef  =  47.64  km, 

Afm  =  8  ■  1010  m4  sec-1, 

Afm  =  3-10-8  sec-1, 

Di  —  600  m, 
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U  =  1  m  sec  . 

Note  that  all  the  dimensional  parameters  are  as  in  McCalpin  and  Haidvogel  [25], 
except  for  the  value  of  the  interfacial  friction  coefficient,  where  they  used  Afm  = 
10“7  sec-1.  With  the  above  values  of  dimensional  parameters,  the  strength  of  the 
total  inflow  into  the  model  domain  is  estimated  at  Tin  —  16  Sv.  After  the  initial  spin 
up,  the  model  was  integrated  for  970  years  to  obtain  all  the  necessary  statistics.  The 
dimensional  time  step  was  140  minutes,  and  the  streamfunction  fields  were  recorded 
every  4.7  days. 

4.4.2  Energy  Time  Series  and  Their  Spectra 

Figure  4-20  shows  the  time  series  of  kinetic,  potential  and  total  energy  which  are 
characterized  by  chaotic  fluctuations.  On  average,  total  energy  is  partitioned  between 
92%  potential  and  8%  kinetic  energy.  Therefore,  total  and  potential  energy  time  series 
are  well  correlated,  and  we  shall  concentrate  on  the  total  energy  behavior  in  the  next 
sections. 

The  kinetic  energy  power  density  spectrum  (Figure  4-21a)  is  red  and  has  dif¬ 
ferent  slopes  for  periods  less  than  4  months  and  periods  between  1  and  7  years.  There 
are  peaks  in  the  spectrum  at  the  time  scales  of  5.1  months  and  7.5  months.  To  en¬ 
sure  their  significance,  we  narrowed  confidence  intervals  by  computing  power  density 
spectrum  for  periods  less  than  5  years  (Figure  4-23b).  The  energy  associated  with 
the  5-month  period  is  at  least  twice  as  high  as  the  energy  of  the  7-month  period 
(Figure  4-23d).  The  first  peak  is  associated  with  the  interaction  of  large  amplitude, 
westward  propagating  meanders  with  the  western  boundary  during  the  low  energy 
period,  which  leads  to  large  meridional  excursions  in  the  separation  point  and  to  for¬ 
mation  of  rings  on  the  wall.  Figure  4-21b  presents  power  density  spectrum  multiplied 
by  the  frequency  which  reveals  that  most  of  kinetic  energy  lies  in  the  band  between  1 
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Figure  4-21:  (a):  Kinetic  energy  power  density  spectrum.  Dashed  line  denotes  con¬ 
fidence  interval  at  95%.  (b):  Frequency  times  power  density  spectrum  for  kinetic 
energy.  Parameters:  Tin  =  16  Sv,  Afm  -  8-1010  m4  sec-1  and  Afm  =  3.0- 10“8  sec-1. 
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Figure  4-22:  (a):  Potential  energy  power  density  spectrum.  Dashed  line  denotes 
confidence  interval  at  95%.  (b):  Frequency  times  power  density  spectrum  for  potential 
energy.  Parameters:  Tin  =  16  Sv,  Afm  =  8-1010  m4  sec-1  and  Afm  —  3.0 -10~8  sec-1. 
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Figure  4-23:  Potential  (a)  and  kinetic  (b)  energy  power  density  spectrum  for  peri¬ 
ods  less  than  5  years.  Corresponding  frequency  times  power  density  spectrum  for 
potential  (c)  and  kinetic  (d)  energy.  Dashed  line  denotes  confidence  interval  at  95%. 
Parameters:  Tin  =  16  Sv,  Afm  =  8  ■  IO10  m4  sec-1  and  Aiim  =  3.0  •  10~8  sec-1. 


and  10  years.  For  periods  higher  than  7  years,  the  kinetic  energy  spectrum  is  roughly 
flat. 

Figure  4-22  shows  that  most  of  the  potential  energy  is  on  the  decadal  time 
scales,  with  a  significant  level  on  interannual  scales.  Similar  to  the  kinetic  energy 
spectrum,  it  has  peaks  at  5.1  and  7.5  months  (Figure  4-23a),  however,  the  energy 
associated  with  them  is  quite  low  (Figure  4-23c),  which  might  be  the  signature  of  the 
ring  life  cycle.  For  comparison,  a  peak  in  the  potential  energy  spectrum  (not  in  the 
kinetic  energy)  at  periods  between  7  and  8  months  was  also  observed  in  the  wind- 
driven  reduced-gravity  model  of  McCalpin  and  Haidvogel  [25]  and  was  attributed  to 
the  preferred  time  scale  of  eddy-mean  flow  interaction. 

4.4.3  Energy  Histograms 

The  histogram  of  nondimensional  total  energy,  TE,  (Figure  4-24)  has  a  bimodal  dis¬ 
tribution  which  reveals  two  preferred  energy  regimes  of  the  model:  state  L  with 
0.062  <  TE  <  0.064  and  state  H  with  0.066  <  TE  <  0.072.  The  transition 
states,  HL  (high  to  low)  and  LH  (low  to  high),  are  characterized  by  the  energy 
level  0.064  <  TE  <  0.066.  As  we  show  in  the  following  section,  HL  and  LH  transi¬ 
tions  are  dynamically  different.  Based  on  the  histogram,  the  probability  of  the  model 
staying  in  a  particular  energy  regime  is  estimated  as  0.79  for  state  H,  0.16  for  state 
L  and  0.05  for  the  transitions. 

For  comparison,  in  the  wind-driven  reduced-gravity  model  of  McCalpin  and 
Haidvogel  [25],  the  total  energy  histogram  did  not  have  a  bimodal  distribution.  In¬ 
stead,  it  displayed  three  energy  peaks  associated  with  three  different  dynamical  regimes 
There  was  significant  overlap  between  the  different  states,  which  indicated  that  the 
model  was  spending  as  much  time  in  transition  as  in  any  given  state. 

The  distribution  of  kinetic  energy  is  not  normal.  MATLAB  QQPLOT  rou¬ 
tine  applied  to  the  kinetic  energy  time  series  shows  the  quantiles  of  the  normalized 
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Figure  4-24:  Histograms  of  kinetic  (a),  potential  (b)  and  total  (c)  energy  for  the 
reference  run.  Parameters:  Tin  =  16  Sv,  A6  =  8  •  1010  m4  sec-1  and  A2  =  3.0  • 
10-8  sec-1. 
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Figure  4-25:  QQPLOT  of  total  kinetic  energy.  Dashed  line  is  error  bars  for  normal 
distribution.  Solid  line  denotes  quantiles  of  the  normalized  kinetic  energy  time  series 
vs  quantiles  of  the  generated  time  series  with  normal  distribution.  Parameters:  T{n  = 
16  Sv,  A§  =  8  •  1010  m4  sec-1  and  A2  =  3.0  •  10-8  sec-1. 

kinetic  energy  time  series,  versus  the  quantiles  of  the  generated  random 

time  series  with  normal  distribution.  Here,  Kj  —  ( KEj  —  KE)/std(KE),  where 
KE  —  jjE^=lKEj  and  std (KE)  are  the  mean  and  the  standard  deviation  of  the  ki¬ 
netic  energy  time  series,  {KEf=1}  .  QQPLOT  allows  to  determine  whether  the  time 
series  come  from  the  same  distribution.  Figure  4-25  indicates  that  the  tails  of  the 
kinetic  energy  time  series  do  not  fit  within  the  error  bars  calculated  for  normal  dis¬ 
tribution.  Therefore,  kinetic  energy  time  series  is  not  normal  due  to  the  existence 
of  two  different  dynamical  regimes.  However,  because  of  a  significant  overlap,  the 
distribution  does  not  appear  to  be  bimodal. 

As  we  describe  later,  in  the  present  model,  energy  state  L  is  related  to  a  dy¬ 
namical  regime  of  an  eastward  flowing  jet  with  strong  meridional  fluctuations  in  the 
position  of  the  separation  point  and  violent  large-amplitude,  westward  propagating 
meanders.  The  state  does  not  have  recirculating  gyres,  and  its  eddy  field  is  char¬ 
acterized  by  strong  rings  generated  only  on  the  western  wall  and  not  in  the  model 
interior.  State  H,  on  the  other  hand,  possesses  well-developed  recirculations.  How¬ 
ever,  the  longitudinal  size  and  the  symmetry  of  the  recirculations  suffer  significant 
changes  in  time,  thus  giving  rise  to  strong  low  frequency  variability  within  energy 
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state  H.  The  eddy  field  is  characterized  by  periods  of  strong  and  weak  ring  forma¬ 
tion  in  the  region  east  of  the  recirculation  gyres.  An  attempt  to  determine  clear 
energetic  boundaries  between  solutions  with  large,  antisymmetric  recirculations  and 
small,  asymmetric  ones  failed.  Therefore,  the  energy  histogram  provides  energetic 
boundaries  only  between  a  state  with  no  recirculations  and  a  state  with  spatially 
varying  recirculations. 

4.4.4  Time  Mean  Flow 

Although  the  total  energy  has  a  bimodal  distribution,  the  time  averaged  streamfunc- 
tion  and  potential  vorticity  fields  are  still  useful  in  describing  the  system.  Averaged 
over  the  whole  time  series,  the  time  mean  streamfunction  field  (Figure  4-26a)  is  char¬ 
acterized  by  two  antisymmetric  recirculating  gyres  pressed  against  the  eastward  jet 
and  the  western  boundary  currents.  The  meridional  and  longitudinal  sizes  of  the 
gyres  are  200  and  1000  km  respectively.  East  of  the  recirculations,  the  jet  widens  up 
to  400  km. 

The  corresponding  potential  vorticity  field  (Figure  4-26b)  shows  two  tongues 
of  anomalously  high  and  low  potential  vorticity  advected  by  the  western  boundary 
current  towards  the  middle  of  the  domain  from  the  north  and  south.  The  recirculation 
regions  are  marked  by  closed  potential  vorticity  contours.  The  value  of  homogenized 
potential  vorticity  in  the  northern  recirculation  indicates  that  water  is  brought  from 
the  north,  entrained  into  the  jet  and  then  mixed  north  of  the  jet.  The  opposite  is 
true  of  the  southern  recirculation.  The  potential  vorticity  front  associated  with  the 
eastward  jet  is  sharp  near  the  separation  point  and  becomes  wider  east  of  500  km, 
which  might  be  the  sign  of  enhanced  cross-front  mixing  in  this  area. 
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Figure  4-26:  Contours  of  the  time  mean  streamfunction  field  (a)  and  potential  vortic- 
ity  field  (b).  Dashed  line  denotes  negative  values,  solid  line  denotes  positive  values. 
Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1  and  Afm  =  3.0  •  10-8  sec-1. 
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4.4.5  Time  Mean  Eddy  Variability 

The  spatial  structure  of  the  time  mean  eddy  variances  (Figures  4-27a-d)  indicates  that 
the  variability  is  concentrated  in  the  area  of  the  eastward  jet  and  its  recirculations  and 
in  the  western  boundary  region  close  to  the  separation  point.  Within  the  eastward 
jet,  the  variability  reaches  its  maximum  around  800-1000  km  and  decays  further 
eastward  leveling  to  zero  beyond  2000  km.  Another  maximum  is  aligned  with  the 
western  boundary  current  between  900  and  1500  km.  As  will  be  discussed  later,  the 
conical  shape  of  the  variability  region  is  associated  with  the  barotropic  instabilities 
in  the  form  of  waves,  rings  and  meanders  and  their  interactions  with  the  western 
boundary  current,  changes  in  the  spatial  structure  and  the  size  of  the  recirculation 
gyres  and  meridional  shifts  in  the  position  of  the  separation  point  of  the  eastward  jet. 

4.4.6  Statistical  Differences  in  the  Eddy  and  Time  Mean 
Flow  Energy 

Figures  4-28a,b  present  a  segment  of  the  total  energy  and  the  corresponding  eddy 
energy  time  series.  Visually,  the  low  total  energy  periods  are  correlated  with  the 
periods  of  maximum  eddy  energy.  To  quantify  it,  we  computed  a  two-dimensional 
histogram  (Figure  4-29)  showing  the  number  of  times  the  model  had  a  particular 
total  energy/eddy  energy  pair.  When  interpreting  the  results,  one  must  remember 
that  total  energy  is  the  sum  of  mean  flow  energy,  eddy  energy  and  temporary  posi¬ 
tive  or  negative  correlations  between  the  mean  and  fluctuating  fields.  The  histogram 
shows  two  separate  clusters  of  the  eddy  energy  values.  The  first  one  is  associated 
with  the  high  total  energy  state,  H,  and  has  the  following  mean  and  standard  devia¬ 
tion:  (m£r,<T/f)  =  (0.0105,0.037).  The  second  cluster  corresponds  to  L-state  and  has 
K,^)  =  (0.0163,0.041)  as  mean  and  standard  deviation.  The  above  analysis  sug¬ 
gests  that  the  high  (low)  total  energy  regimes  are  associated  with  the  low  (high)  eddy 
energy.  The  small  overlap  between  the  two  eddy  states  indicates  that  in  addition  to 
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Figure  4-28:  A  segment  of  the  time  series  of  nondimensional  total  energy  (a)  and 
of  corresponding  nondimensional  total  eddy  energy  (b).  Parameters:  Tin  =  16  Sv, 
Afm  =  8  •  1010  m4  sec”1  and  Afm  =  3.0  •  lO^sec"1. 
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Figure  4-29:  Two-dimensional  histogram  of  a  segment  of  the  total  energy  and  the 
total  eddy  energy  from  Figure  4-28.  Parameters:  T{n  =  16  Sv,  Af*m  =  8  •  1010  m4 
sec-1  and  Af™  =  3.0  - 10-8  sec-1. 

changes  in  the  eddy  intensity,  there  might  be  some  structural  differences  in  the  eddy 
field. 


4.4.7  Enstrophy 


Figure  4-30a  presents  the  histogram  of  total  nondimensional  enstrophy  for  the  ref¬ 
erence  run.  Clearly,  similar  to  total  energy,  it  has  a  bimodal  distribution  with  well- 
separated  peaks  associated  with  two  dynamical  regimes.  Hence,  in  addition  to  total 
energy,  total  enstrophy  could  be  another  indicator  of  different  dynamical  regimes. 
The  distribution  of  perturbation  enstrophy  (Figure  4-30b),  although  also  bimodal, 
has  a  significant  overlap. 

Total  enstrophy  length  scale  can  be  defined  as: 


(4.11) 


where  L  =  2400  km  is  model’s  horizontal  length  scale,  E  is  domain-integrated  total 

2 

nondimensional  energy  and  e  =  J^n  XT'  \  ls  total  nondimensional  enstrophy. 


147 


g  600 


8 400 
o 

£  300 
O  200 

z 

UJ 

o100 

UJ 

£  0 
ie 

g  600 
z 

UJ 

tr  500 
a: 

§400 
o 

g  300 

O  200 
z 

UJ 

8 100 
UJ 

£  o 

i 

Figure  4-30:  (a):  Histogram  of  total  nondimensional  enstrophy.  (b)  Histogram  of 
perturbation  nondimensional  enstrophy.  Parameters:  Tin  —  16  Sv,  Afm  =  8  •  1010  m4 
sec-1  and  Afm  =  3.0  •  10-8  sec-1. 

Similarly,  eddy  enstrophy  length  scale,  L'enst,  is  defined  in  (4.11),  where  E  is  now  the 
domain-integrated  perturbation  energy  and  e  =  ^  dydx  is  eddy  enstrophy. 

Figure  4-31  shows  two-dimensional  histograms  of  total  energy  vs  Lenst  and  of 
eddy  energy  vs  L'enst.  Clearly,  there  is  no  significant  change  in  the  enstrophy  length 
scale  between  the  high  and  low  energy  regimes.  Therefore,  regime  variability  is  not 
connected  to  scale  variability. 

4.4.8  Mean  Fields  Associated  with  Different  Energy  States 

The  averaging  of  the  fields  within  several  energy  bands  shows  qualitative  differences 
in  the  structure  of  the  flow  and  helps  to  understand  the  nature  of  the  multiple  dy¬ 
namical  states  of  the  model.  Figures  4-32  and  4-33  present  ensemble  averages  of  the 
streamfunction  and  potential  vorticity  over  four  energy  bands:  low  energy:  L=[0.062 
0.064],  transition  and  a  low  range  of  the  high  energy  state:  HL=[0.064,  0.0685]  ,  a 
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Figure  4-31:  (a):  Two-dimensional  histogram  of  nondimensional  total  energy,  E,  vs 
dimensional  total  enstrophy  scale,  Lenst  =  EJ  2',  where  q  is  potential  vortic- 
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ity.  (b):  Two-dimensional  histogram  of  nondimensional  perturbation  energy,  E,  vs 
dimensional  perturbation  enstrophy  scale,  Lensi  =  LyJ — ^ ,  where  q  is  pertur¬ 
bation  potential  vorticity.  Parameters:  Tjn  =  16  Sv,  =  8  •  1010  m4  sec-1  and 
Afm  =  3.0  ■  10“8  sec"1,  L  =  2400  km. 
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Figure  4-32:  Ensemble  average  of  the  streamfunction  fields  with  the  total  energy  levels 
for  HH-state  (TE  >  0.0705)  (a),  HM-state  (0.0685  <  TE  <  0.0705)  (b),  HL-state 
(0.064  <  TE  <  0.0685)  (c)  and  L-state  (TE  <  0.064)  (d).  Parameters:  T*n  =  16  Sv, 
Afm  =  8  •  1010  m4  sec-1  and  Af171  =  3.0  •  10-8  sec-1. 
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Figure  4-33:  Ensemble  average  of  potential  vorticity  fields  with  the  total  energy  levels 
for  HH-state  (TE  >  0.0705)  (a),  HM-state  (0.0685  <TE  <  0.0705)  (b),  HL-state 
(0.064  <  TE  <  0.0685)  (c)  and  L-state  (TE  <  0.064)  (d).  Parameters:  Tin  =  16  Sv, 
Afm  =  8  •  1010  m4  sec-1  and  Afm  =  3.0  •  10-8  sec-1. 
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medium  range  of  the  high  energy  state:  HM=[0.0685,  0.0705]  and  a  high  range  of 
the  high  energy  state:  HH=[0.0705,  0.072]. 

The  structures  of  the  streamfunction  fields  are  similar  for  the  last  three  energy 
ranges:  all  of  them  are  characterized  by  the  presence  of  symmetric  recirculations. 
However,  the  longitudinal  size  of  the  gyres  decreases  from  1500  km  for  HH-state 
to  800  km  for  HL-state.  In  HM-  and  HL-  states,  the  width  of  the  eastward  jet 
increases  from  200  km  just  east  of  the  recirculations  to  400  km  east  of  2000  km.  The 
sudden  change  in  the  jet’s  width  is  associated  with  the  presence  of  larger  meanders 
and  the  formation  of  rings  east  of  the  recirculation  gyres.  The  corresponding  potential 
vorticity  fields  differ  not  only  by  the  shape  of  the  recirculation  regions,  but  also  by  the 
values  of  the  mixed  q  :  HH  has  the  highest  values  of  anomalous  potential  vorticity. 
Another  obvious  difference  is  the  detachment  of  the  mixed  q— region  from  higher  and 
lower  latitudes  in  HL-state,  which  allows  for  larger  excursion  of  potential  vorticity 
contours  indicating  a  more  unstable  flow  structure. 

In  the  L-energy  range,  the  time  mean  field  is  a  combination  of  equally  possible 
antisymmetric  solutions.  One  is  characterized  by  the  overshooting  southward-flowing 
western  boundary  current  separating  from  the  wall  approximately  300  km  south  of 
the  middle  of  the  domain.  The  other  one  has  a  northward-flowing  western  boundary 
current  separating  300  km  north  of  the  middle  of  the  domain.  Such  solutions  explain 
the  presence  in  the  time  mean  of  an  anticyclonic/cyclonic  vortex  north/south  of  the 
middle  of  the  domain.  In  the  time  mean,  the  recirculations  appear  to  be  pushed 
further  eastward,  and  their  zonal  extent  is  only  500  km.  They  exist  mostly  in  the 
time  mean  sense;  the  instantaneous  fields  have  individual  vortices  confined  between 
large  amplitude,  westward  propagating  meanders.  The  corresponding  time  mean 
potential  vorticity  field  appears  to  have  a  500  km  by  500  km  region  of  “homogenized” 
potential  vorticity  in  the  western  part  of  the  domain.  This  is  misleading,  since  the 
time  mean  field  is  constructed  by  averaging  of  two  opposite  sign  potential  vorticity 
fields  due  to  overshooting  events.  Nevertheless,  this  is  a  region  of  vigorous  mixing 
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across  the  front,  since  large  amplitude,  westward  propagating  meanders  interacting 
with  the  western  wall  form  strong  rings,  that  propagate  away  from  the  middle  of  the 
domain. 

4.4.9  Typical  Energy  Cycle 

Eddy  Energy  -  Total  Energy  Partition 

To  understand  better  the  nature  of  multiple  dynamical  regimes,  we  choose  a  typical 
cycle  (Figure  4-34a)  when  the  total  energy  rapidly  decays  (HL-transition)  and,  after 
fluctuating  around  minimum  values  for  10-50  years,  quickly  builds  up  again  reaching 
its  maximum  (LH-transition) .  Typical  transitions  take  between  5  and  10  years  and 
are  associated  with  significant  changes  in  the  spatial  structure  of  the  mean  and  eddy 
fields.  Such  changes  are  dynamically  asymmetric  for  HL-  and  LH-transitions.  The 
ratio  of  eddy  to  total  energy  (Figure  4-34b)  indicates  that  in  the  former  case,  the  eddy 
field  part  increases  from  10%  to  30%.  During  the  LH-transition,  the  eddy  energy  level 
reaches  its  minimum  of  7%  and  averages  at  15%  in  H-state.  The  amplitude  of  the 
fluctuations  in  the  eddy-mean  flow  partition  is  the  highest  during  L-state. 

Flow  Evolution  During  a  Typical  Energy  Cycle 

In  this  section,  we  describe  the  evolution  of  the  streamfunction  and  potential  vortic- 
ity  fields  during  a  typical  energy  cycle  depicted  in  Figure  4-34a.  At  t  =  1.3  years 
(denoted  by  number  1  in  Figure  4-34a),  when  the  total  energy  is  in  H-state,  the 
flow  (Figures  4-35a,b)  has  well-developed  recirculations  with  a  longitudinal  extent 
of  1500  km.  Small  fluctuations  in  the  position  of  the  separation  point  are  due  to 
weak  eastward  propagating  meanders.  East  of  the  recirculation  region,  where  the  jet 
is  wider  and  therefore  more  stable,  meanders  are  characterized  by  westward  propa¬ 
gation  and  larger  amplitudes.  The  collision  of  eastward  and  westward  propagating 
meanders  leads  to  formation  of  rings  in  the  region  east  of  the  recirculations.  The 
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Figure  4-34:  (a):  Segment  of  the  total  energy  time  series  corresponding  to  a  typical 
low  frequency  cycle,  (b):  The  corresponding  ratio  of  eddy  to  total  energy  times  100%. 
Parameters:  Tin  =  16  Sv,  A  =  8  •  1010  m4  sec"1  and  A  =  3.0  •  10-8  sec-1. 


meander  amplitudes  define  the  strength  and  the  size  of  created  vortices  and  thus  the 
degree  of  potential  vorticity  mixing  across  the  front.  Once  the  ring  detaches  from  the 
front,  it  travels  westward  in  the  recirculating  gyres  mixing  the  opposite  sign  vorticity 
into  the  surroundings  and  homogenizing  the  surrounding  gradient.  Upon  reaching 
the  western  wall,  at  this  stage  of  the  flow  development,  rings  are  not  strong  enough 
to  overcome  the  advective  influence  of  the  western  boundary  current.  Thus,  they  join 
the  eastward  jet  and  further  perturb  it. 

At  t  =  5.2  years  (denoted  by  number  2  in  Figure  4-34a),  the  total  energy  level 
is  still  high,  and  the  eddy  energy  partition  increases  from  10%  to  20%.  The  growth 
of  jet  instabilities  and  subsequent  ring  formation  lead  to  a  decrease  to  800  km  in 
the  longitudinal  extent  of  the  recirculations  (Figures  4-35c,d).  Although  the  jet  front 
is  still  symmetric  in  the  recirculation  region,  the  gyres  lose  their  symmetry.  One  of 
them  starts  wrapping  around  another  one,  thus  giving  birth  to  stronger  and  larger 
rings  and  therefore  providing  enhanced  potential  vorticity  mixing  across  the  front. 
Rings  in  turn  tend  to  further  decrease  the  size  and  the  strength  of  the  recirculations. 

Time  t  =  10.4  years  (denoted  by  number  3  in  Figure  4-34a)  marks  the  begin¬ 
ning  of  a  sharp  decline  in  the  total  energy  level  and  corresponds  to  a  further  increase 
in  the  eddy  energy  partition.  Although  the  jet  still  separates  in  the  middle  of  the 
domain,  its  front  exhibits  strong  northward  and  southward  oscillations  east  of  500 
km.  Figures  4-35e,f  show  one  of  the  states,  with  a  smaller  northern  recirculation 
and  a  northward  front  excursion.  The  latter  results  in  a  production  of  stronger  and 
larger  rings  on  the  south  and  weaker  and  smaller  rings  on  the  north.  Therefore,  upon 
reaching  the  western  wall,  the  rings  have  different  trajectories;  strong  cyclonic  rings 
are  propelled  southward,  thereby  decreasing  the  inflow  from  the  south,  whereas  weak 
anticyclonic  rings  join  the  eastward  jet  without  affecting  the  inflow  from  the  north 
to  any  great  extent.  These  processes  create  conditions  favorable  for  the  southward 
shift  in  the  position  of  the  separation  point.  The  opposite  is  true  when  the  front  has 
a  southward  excursion. 
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Figure  4-36  presents  instantaneous  streamfunction  and  potential  vorticity  fields 
at  times  t  =  13,  14.3  and  15.6  years  (denoted  by  numbers  4,  5  and  6  in  Figure  4-34a) 
which  correspond  to  a  30%  maxima  in  the  eddy  energy  partition.  A  sharp  decay  in 
the  total  energy  level  is  due  to  the  potential  energy  decline  associated  with  an  almost 
complete  disappearance  of  the  recirculations.  Instead  of  the  well-developed  gyres, 
there  are  individual  vortices  corresponding  to  meander  peaks  and  troughs  east  of  500 
km.  Such  structures  provide  efficient  mixing  on  the  flanks  of  an  eastward  jet.  The 
separation  point  shifts  300  km  south  and  north  of  the  middle  of  the  domain  and, 
along  with  large  amplitude  meanders  provides  stronger  fluctuations  in  the  domain- 
integrated  potential  vorticity  field. 

During  this  stage  of  the  flow  development,  there  is  no  ring  formation  in  the 
interior  of  the  domain.  Instead,  strong  rings  are  constantly  formed  on  the  western 
boundary  in  the  following  manner.  Large  amplitude  meanders  propagate  westward 
similar  to  long  Rossby  waves.  If  a  negative  potential  vorticity  anomaly  associated 
with  a  meander  peak  reaches  the  western  wall,  it  starts  moving  northward  due  to  the 
influence  of  an  image  anomaly  arising  from  a  no-normal-flow  boundary  condition.  The 
separation  point  is  therefore  shifting  northward.  Meanwhile,  an  upstream  positive 
potential  vorticity  anomaly  corresponding  to  a  meander  trough  continues  to  propagate 
westward.  Upon  hitting  the  western  wall,  it  starts  moving  southward.  Two  anomalies 
propelling  in  the  opposite  directions  cause  the  disruption  in  the  potential  vorticity 
front.  As  a  result,  a  strong,  northward  propagating  anticyclonic  ring  is  formed  from 
the  negative  potential  vorticity  anomaly,  and  the  separation  point  jumps  southward. 

The  process  is  repeated  for  the  positive  potential  vorticity  anomaly,  which 
leads  to  further  southward  excursion  of  the  separation  point  and  to  the  formation  of  a 
strong,  southward  propagating  cyclonic  ring.  The  dynamics  of  such  strong  boundary- 
formed  rings  tends  to  restore  the  symmetry  of  the  inflow  conditions.  However,  the 
interaction  of  the  large  amplitude,  westward  propagating  meanders  leads  to  the  per¬ 
sistence  of  a  low  energy  state. 
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The  model  stays  in  this  dynamical  regime  until  the  instabilities  start  decaying, 
perhaps  due  to  frictional  effects,  which  results  in  the  decrease  of  meander  amplitudes 
and  in  the  amplitude  of  the  separation  point  excursions.  The  strength  of  the  boundary 
formed  rings  also  decreases,  which  helps  to  restore  the  inflow  symmetry.  At  times 
t  =  18.2  and  20.8  years  (Figures  4-37a-d),  which  corresponds  to  the  middle  and 
to  the  end  of  state  L  (denoted  by  numbers  7  and  8  in  Figure  4-34a),  the  shifts 
of  the  separation  point  become  so  small  that  weak  recirculations  develop  near  the 
western  boundary.  At  this  point,  rings  form  just  east  of  the  recirculations,  and  their 
detachment  tends  to  flatten  the  front. 

The  transition  to  the  high  energy  state  is  illustrated  by  Figures  4-37e-f  and  4- 
38a-d,  which  are  snapshots  at  t  =  21,  23.5,  26.1  and  28.7  years  (denoted  by  numbers 
9-12  in  Figure  4-34a).  The  eddy  energy  partition  drops  below  10%  during  this  period. 
The  size  of  the  recirculations  increases  from  500  to  almost  2000  km,  and  the  ring 
formation  process  slows  down  and,  at  some  point,  halts.  Thus,  the  model  reaches 
the  dynamical  state  associated  with  large  recirculation  gyres  and  very  weak  eddy 
activity.  The  spatial  structure  resembles  the  antisymmetric  steady  solution  described 
in  section  4.2.1. 

However,  the  energy  does  not  stay  in  its  maxima  for  a  long  time.  Instead, 
the  model  trajectory  reaches  the  vicinity  of  an  antisymmetric  steady  solution  and 
then  quickly  escapes  back  to  the  regime  where  nonlinear  eddy-eddy  interactions  are 
important.  This  might  happen  because  the  model  has  no  mechanism  that  would 
restore  the  symmetry  of  the  inflow  conditions  near  the  separation  point  and  therefore 
guide  the  model  trajectory  towards  the  fixed  point.  In  the  wind-driven  models  with 
nearly  antisymmetric  wind-stress  curl,  the  correlation  between  the  wind-stress  curl 
and  the  antisymmetric  steady  solution  might  be  important  in  helping  the  model 
trajectory  to  reach  a  closer  neighborhood  of  the  fixed  point.  In  agreement  with  the 
linear  stability  analysis  of  Primeau  [29] ,  it  takes  more  than  20  years  for  the  instabilities 
to  grow  and  to  move  the  wind-driven  model  from  a  high  to  a  medium  energy  state. 
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Figure  4-35:  Contours  of  instantaneous  streamfunction  field  and  potential  vorticity 
at  time  t  =  1.3  years  (a,b),  t  =  5.2  years  (c,d)  and  t  =  10.4  years  (e,f).  Time, 
t,  is  as  in  Figure  4-34.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec"1  and 
Afm  =  3.0  •  10-8  sec-1.  The  snapshots  are  presented  for  the  western  part  of  the 
model  domain. 
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Figure  4-36:  Contours  of  instantaneous  streamfunction  field  and  potential  vorticity 
at  time  t  =  13.05  years  (a,b),  t  =  14.35  years  (c,d)  and  t  =  15.6  years  (e,f).  Time, 
t,  is  as  in  Figure  4-34.  Parameters:  Tin  =  16  Sv,  Afm  —  8  •  1010  m4  sec-1  and 
A f"1  =  3.0  •  10-8  sec-1.  The  snapshots  are  presented  for  the  western  part  of  the 
model  domain. 


159 


(a)  (b) 


(c) 


(d) 


(e) 


(f) 


Figure  4-37:  Contours  of  instantaneous  streamfunction  field  and  potential  vorticity 
at  time  t  =  18.2  years  (a,b),  t  =  20.8  years  (c,d)  and  t  =  23.5  years  (e,f).  Time, 
t,  is  as  in  Figure  4-34.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1  and 
Af171  =  3.0  •  10-8  sec-1.  The  snapshots  are  presented  for  the  western  part  of  the 
model  domain. 
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Figure  4-38:  Contours  of  instantaneous  streamfunction  field  and  potential  vorticity 
at  time  t  =  26.1  years  (a,b),  t  =  28.7  years  (c,d)  and  t  =  32.6  years  (e,f).  Time, 
t,  is  as  in  Figure  4-34.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1  and 
Afm  =  3.0  •  10-8  sec-1.  The  snapshots  are  presented  for  the  western  part  of  the 
model  domain. 
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Figure  4-39:  Normalized  histograms  of  domain-integrated  streamfunction  (a,  solid) 
and  potential  vorticity  (b,  solid)  for  the  reference  run.  Dashed  line  denotes  histograms 
for  the  first  half  of  the  record,  dashed-dotted  line  denotes  histograms  for  the  second 
half  of  the  record.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec"1  and  Afm  = 
3.0  •  10“8  sec"1. 

4.4.10  Mass  and  Potential  Vorticity  Balances 

Before  we  describe  the  mass  and  potential  vorticity  balances,  let  us  recall  that  due  to  a 
hydrostatic  balance  and  quasi-geostrophic  assumptions,  mass  is  given  by  the  domain- 
integrated  streamfunction  field.  The  present  model  does  not  have  special  numerical 
procedures  which  would  guarantee  conservation  of  mass  and  potential  vorticity  at  each 
time  step.  However,  fixed  boundary  conditions,  providing  the  balance  of  geostrophic 
inflows  and  outflows,  and  the  quasi-geostrophic  equations  of  motion,  allowing  only  for 
an  order  of  Rossby  number  disturbances  in  the  pressure  surfaces,  permit  only  small 
deviations  in  the  domain-integrated  streamfunction  from  its  original  zero  value. 

Figure  4-39  shows  normalized  histograms  of  the  domain-integrated  streamfunc¬ 
tion  and  potential  vorticity  fields.  Dashed  and  dashed-dotted  lines  present  similar 
histograms  based  on  the  first  and  the  second  halves  of  the  time  series  and  indicate 
that  the  shape  of  the  histograms  does  not  change  qualitatively.  Therefore,  the  length 


162 


of  the  time  series  is  pretty  robust.  A  straight  line  fitted  to  the  time  series  gives  slope 
coefficients  of  6.1e-10  and  l.le-6  for  the  integrals  of  the  streamfunction  and  potential 
vorticity  fields  respectively.  Although  the  realization  is  finite,  the  small  values  of  the 
slope  coefficients  indicate  that  integrated  mass  and  potential  vorticity  do  not  drift 
in  time,  even  though  the  model  has  no  explicit  conservation  scheme  for  any  of  those 
fields. 

QQPLOTs  based  on  the  time  series  of  domain-integrated  streamfunction  and 
potential  vorticity  are  presented  in  Figures  4-40a,c.  Clearly,  there  are  large  deviations 
from  a  linear  function,  which  indicates  that  the  corresponding  distributions  are  not 
normal.  The  latter  is  due  to  the  presence  of  the  low  energy  state.  Figure  4-41  shows 
histograms  based  on  the  segments  of  the  original  time  series  of  fx  f  ip  dxdy  and 
fx  f  q  dxdy  that  were  constructed  when  total  energy  was  either  larger  than  0.067 
(a  high  energy  state,  solid  line)  or  smaller  than  0.064  (a  low  energy  state,  dashed 
line).  The  histograms  of  ip—  and  q— integrals  have  normal  distributions  during  the 
high  energy  period,  since  the  corresponding  QQPLOTs  are  linear  (Figures  4-40b,d). 
As  we  saw  earlier,  during  the  high  energy  state,  the  flow  structure  is  characterized  by 
the  well-developed  recirculations,  and  the  time  mean  streamfunction  and  potential 
vorticity  fields  are  well  defined. 

However,  during  the  low  energy  state,  the  model  displays  strong  fluctuations 
associated  with  large  meridional  excursions  in  the  separation  point.  Clearly,  when 
there  is  significant  northward  (southward)  excursion,  the  integrated  potential  vorticity 
field  must  be  negative  (positive).  The  distributions  of  ip—  and  q— integrals  during  the 
low  energy  period  are  not  normal,  since  there  is  no  well  defined  time  mean  and  since 
there  is  no  preferred  flow  structure.  Therefore,  the  presence  of  the  low  energy  state 
is  responsible  non-normal  distribution  of  the  mass  and  integrated  potential  vorticity 
time  series. 
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Figure  4-40:  (a):  QQPLOT  of  the  quantiles  of  the  time  series  of  fxfyip  dxdy  vs 
the  quantiles  of  the  randomly  generated  time  series  with  normal  distribution,  (b): 
QQPLOT  of  the  time  series  of  Jx  fy  ip  dxdy  with  the  corresponding  total  energy  level 
greater  than  0.067.  (c):  QQPLOT  of  the  total  time  series  of  the  domain-integrated 
potential  vorticity.  (d):  QQPLOT  of  the  domain-integrated  potential  vorticity  time 
series  with  the  corresponding  total  energy  level  greater  than  0.067.  Dashed  line 
denotes  error  bars  for  normal  distribution.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010 
m4  sec-1  and  Afm  =  3.0  •  10-8  sec-1. 
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Figure  4-41:  Normalized  histograms  of  domain-integrated  streamfunction  (a)  and 
potential  vorticity  (b)  for  the  reference  run.  Dashed  line  denotes  histograms  based 
on  the  segment  of  the  total  time  series  when  total  energy  level  was  larger  than  0.067. 
Dashed  line  corresponds  to  the  segments  when  the  total  energy  level  was  smaller  than 
0.064.  Parameters:  Tin  =  16  Sv,  Afm  =  8  •  1010  m4  sec-1  and  Af™  =  3.0  •  10-8  sec-1. 
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4.4.11  Transitions  Between  High  and  Low  Energy  States 

To  understand  the  dynamics  of  transitions  between  the  high  and  low  energy  states, 
the  “ensemble  average”,  or  composite,  time  series  of  the  streamfunction  field  was 
constructed  in  the  following  way.  First,  we  identified  all  the  HL  (LH)  events.  Second, 
we  chose  subjectively  a  “mean”  energy  level  for  all  transitions  as  that  at  the  middle 
of  the  transition.  Third,  we  extracted  2.5  years  of  ^-field  before  and  after  the  model 
reached  this  energy  value.  The  latter  permitted  time  series  of  equal  lengths  for  HL 
(LH)  events: 

{^fc(h),  ^k(h),  •••,  ipk(tm)}, 

where  m  is  the  length  of  each  time  series  and  k  =  1,  •  •  • ,  K  is  the  index  of  an  individual 
event.  Third,  a  “mean”  HL  (LH)  transition  time  series  was  constructed  as 

•••  i/>(tm)},  (4.12) 

where  ^(ti)  =  £E*Li  ^k{U). 

Figures  4-42a-b  show  the  “time  mean”  streamfunction  and  potential  vorticity 
fields  calculated  from  the  time-series  defined  by  (4.12)  for  the  HL-transitions.  One 
should  remember  that  since  transitions  are  fundamentally  non-steady,  the  term  “time 
mean”  should  be  used  with  caution.  The  spatial  structures  of  the  fields  are  similar  to 
those  of  the  low  energy  periods;  the  western  boundary  currents  appear  separating  300 
km  north  and  south  of  the  middle  of  the  domain,  and  weak  recirculations  are  shifted 
eastward.  The  above  mean  state  does  not  reflect  the  true  structure  of  the  flow,  when 
the  interaction  of  large-amplitude,  westward  propagating  meanders  with  the  western 
wall  produces  large  meridional  excursions  in  the  separation  point.  The  structure  of 
the  time  mean  perturbation  fields  (Figure  4-43)  reveals  that  dominant  eddy  processes 
and  conversions  occur  along  the  western  boundary. 

During  the  LH-transitions,  the  symmetry  of  the  inflow  near  the  collision  point 
is  restored,  and  the  time  mean  streamfunction  and  potential  vorticity  (Figures  4- 
42c-d)  are  characterized  by  the  1000-km  long  recirculating  gyres.  The  size  of  the 
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Figure  4-42:  Contours  of  mean  streamfunction  and  potential  vorticity  for  HL- 
transition  (a,b)  and  LH-transition  (c-d).  Parameters:  Tin  =  16  Sv,  =  8  •  1010 
m4  sec"1  and  Afm  =  3.0  •  IQ-8  sec"1. 
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Figure  4-43:  Time  mean  perturbation  kinetic  energy  density  (a),  perturbation  po¬ 
tential  energy  density  (b),  barotropic  conversion  energy  density  (c)  and  perturba¬ 
tion  energy  dissipation  (d)  during  the  HL-transition.  Parameters:  Tin  =  16  Sv, 
Afm  =  8  •  1010  m4  sec-1  and  Afm  =  3.0  •  10-8  sec-1. 
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Figure  4-44:  Time  mean  perturbation  kinetic  energy  density  (a),  perturbation  po¬ 
tential  energy  density  (b),  barotropic  conversion  energy  density  (c)  and  perturba¬ 
tion  energy  dissipation  (d)  during  the  LH-transition.  Parameters:  Tin  =  16  Sv, 
Afm  =  8  •  1010  m4  sec”1  and  A%im  =  3.0  •  10“8  sec"1. 
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recirculations  increases  from  500  to  1800  km  during  the  transition  events.  The  max¬ 
ima  in  the  eddy  density  fields  (Figure  4-44)  are  aligned  with  the  axis  of  the  jet,  and  the 
conical  shape  of  the  perturbation  energy  density  indicates  the  importance  of  rings  in 
flattening  of  the  potential  vorticity  front.  The  model  does  not  quite  reach  the  steady 
antisymmetric  solution,  since  the  amplitude  of  the  eddy  field  starts  increasing  when 
the  recirculations  become  large. 

Therefore,  the  dynamics  of  transitions  from  the  high  to  low  energy  states  are 
drastically  different  from  those  of  the  low  to  high  energy  transitions.  In  the  first 
case,  the  eddy  field  becomes  very  strong,  giving  rise  to  the  large  amplitude,  westward 
propagating  meanders,  which  produce  large  meridional  excursion  in  the  separation 
point  and  strong  northward  and  southward  propagating  rings.  In  the  second  case,  the 
eddy  field  weakens,  and  the  advection  of  potential  vorticity  by  the  western  boundary 
currents  restores  the  flow  symmetry  by  quickly  inducing  large  symmetric  recircula¬ 
tions.  Figures  4-45  and  4-46  summarize  the  energy  conversion  during  the  transition 
events.  When  reading  the  diagrams,  one  should  remember  that  the  flow  is  not  in 
a  steady  state  during  the  transition  events,  and  therefore  the  energy  terms  do  not 
balance  each  other.  Both  eddy  kinetic  and  eddy  potential  energies  are  of  the  same 
order  of  magnitude  as  the  “time  mean”  flow  kinetic  energy  during  the  HL-transition. 
The  positive  sign  of  the  domain-integrated  barotropic  conversion  rate  indicates  that 
eddy  energy  grows  at  the  expense  of  the  mean  flow  energy.  During  the  LH-transition, 
the  eddies  release  energy  back  to  the  mean  flow,  which  is  associated  with  the  negative 
barotropic  conversion  integral. 

4.5  Discussion  of  Nonsymmetric  Inflow  Conditions 

In  this  study,  we  used  symmetric  inflow  and  outflow  conditions,  which  were  con¬ 
structed  to  simulate  a  boundary  with  the  Sverdrup  interior  forced  by  symmetric  wind 
forcing.  It  would  be  interesting  to  see  whether  low  frequency  variability  and  mul- 
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Figure  4-45:  Energy  diagram  for  synthetic  HL-transition.  Parameters:  Tin  =  16  Si 
Afm  _  g  .  1010  m4  sec-l  an(J  =  3.0  •  10-8  sec-1. 


Figure  4-46:  Energy  diagram  for  synthetic  LH-transition.  Parameters:  Tin  =  16  Si 
Afm  =  8  •  1010  m4  sec-1  and  Afm  =  3.0  •  10-8  sec-1. 
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tiple  dynamical  regimes  are  still  possible  if  the  strength  of  the  inflows  through  the 
meridional  boundaries  were  not  equal.  In  this  case,  the  model  solution  might  depend 
strongly  on  the  streamfunction  and  potential  vorticity  distributions  on  the  eastern 
boundary,  as  explained  below. 

One  way  to  illustrate  possible  dependence  of  the  model  solution  on  the  eastern 
boundary  outflow  when  the  meridional  inflows  are  of  different  strengths  is  to  calculate 
the  total  potential  vorticity  input  through  the  model  domain.  Assume  that  the  flow 
through  the  southern  part  of  the  domain  is  weaker  than  that  through  the  northern  part 
of  the  domain,  i.e.  ijjs  =  -aipn,  where  0  <  a  <  1,  =  ip(x,ys )  and  ^>n  =  ^(z,yn). 

The  potential  vorticity  distribution  on  the  southern  boundary,  qs  =  q(x,ys),  can 
be  then  expressed  in  terms  of  the  potential  vorticity  distribution  on  the  northern 
boundary,  qn  =  q(x,yn),  as 


Qs  =  Pys  +  -  7 V.  =  (a  -  1  )(3yn  ~  aqn. 

The  total  potential  vorticity  input  into  the  model  domain  is  given  by 


<p_q  u  •  n 
J  n 


dfi, 


£\a(a-l)(3yn^  +  (l-^)qn 
dy 


d^n 


dx 


dx  — 


(4.13) 


(4.14) 


[V»  dip, 

L  q^dy’ 


where  qe  =  q(xe,y)  and  =  ip(xe,y)  are  potential  vorticity  and  streamfunction 
on  the  eastern  boundary,  u  =  (u,  v)T  is  a  geostrophic  velocity  vector  on  the  model 
boundary,  fi,  and  n  is  a  vector  normal  to  the  boundary.  In  the  case  of  symmetric 
inflow  conditions,  the  coefficient  a  is  equal  to  one,  and  therefore  the  first  integral  of 
the  right  hand  side  is  zero.  If  ipe  is  an  odd  function,  the  second  integral  on  the  right 
hand  side  is  also  zero.  Therefore,  total  potential  vorticity  input  is  zero  if  the  inflow 
and  outflow  conditions  are  symmetric.  This  case  was  of  our  primary  interest. 

For  a  nonzero  a ,  the  value  of  the  first  integral  is  negative.  Hence,  in  order 
for  the  total  potential  vorticity  flux  to  vanish,  the  outflow  on  the  eastern  wall  has  to 
balance  potential  vorticity  input  through  the  meridional  boundaries.  If  the  outflow  is 
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such  that  the  total  q—  input  is  non  zero,  it  may  lead  to  enhanced  potential  vorticity 
transfer  across  the  middle  of  the  domain,  and  therefore  it  might  significantly  affect 
the  interior  solution. 

Another  reason,  why  the  model  might  display  different  behavior,  is  that  steady 
state  solutions  might  change  significantly  when  the  inflows  are  chosen  to  be  nonsym- 
metric.  In  the  present  study,  we  did  not  attempt  to  find  such  solutions.  However, 
as  we  showed  for  the  symmetric  inflows,  steady  state  solutions  are  important  for  the 
existence  of  multiple  dynamical  regimes.  Therefore,  one  might  expect  that  the  change 
in  the  structure  of  steady  solutions  might  lead  to  the  change  in  the  character  and  the 
occurrence  of  multiple  states. 

For  comparison,  McCalpin  and  Haidvogel  [25]  showed  that  for  strongly  asym¬ 
metric  wind,  the  reduced-gravity  model  was  not  able  to  reach  the  high  energy  state 
and  instead  stayed  in  the  low  energy  regime.  The  absence  of  the  high  energy  state 
may  have  been  due  to  “correlation”  phenomena  ( Scott  and  Straub  [34])  between  the 
global  structures  of  the  wind  forcing  and  the  interior  solution.  In  our  boundary-forced 
model,  such  correlations  are  not  possible.  Hence,  compared  to  the  wind-driven  model, 
the  colliding  jets  model  might  display  stronger  low  frequency  variability  in  the  case 
of  nonsymmetric  inflow  conditions,  if  the  steady  states  are  close  to  each  other  in  the 
phase  space. 

4.6  Summary  and  Discussion 

In  this  Chapter,  we  found  steady  antisymmetric  solutions  of  the  reduced-gravity, 
colliding  jets  model.  For  a  wide  range  of  interfacial  friction  coefficients,  the  steady 
flow  consisted  of  an  eastward  jet  flanked  by  two  recirculating  gyres.  For  a  set  of  chosen 
parameters,  the  zonal  extent  of  recirculations  and  the  total  energy  level  increased  with 
a  decrease  in  the  interfacial  friction  coefficient. 
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The  initial  value  technique  was  used  to  calculate  the  exponentially  growing 
solutions  of  the  antisymmetric  states.  The  amplitude  of  the  unstable  mode  closely 
resembled  that  of  the  unstable  solution  of  the  wind-driven  reduced-gravity  model  of 
Primeau  [29]  and  consisted  of  four  strong  cells  sitting  on  the  axis  of  the  jet  between  the 
recirculating  gyres  and  weaker  cells  in  the  recirculation  regions.  The  e-folding  time 
varied  between  4.8  years  and  417  days.  With  the  growth  in  the  instability  amplitude, 
the  linear  growth  halted  due  to  nonlinear  interactions  of  the  cells.  That  interaction 
led  to  formation  of  rings  in  the  area  where  the  jet  widens. 

Fully  nonlinear  calculations  showed  that  with  an  increase  in  the  value  of  the 
interfacial  friction  coefficient,  the  model  became  weakly  unstable  and  displayed  a 
unique  dynamical  regime  corresponding  to  an  eastward  jet  with  two  recirculating 
gyres.  For  near  critical  values  of  the  friction  coefficient,  the  eddy  field  was  dominated 
by  the  eastward  and  westward  propagating  meanders  of  small  amplitudes.  Rings 
were  weak  and  therefore  dissipated  during  their  propagation  through  the  recirculation 
regions.  The  energy  level  and  the  flow  structure  were  close  to  those  of  the  steady 
antisymmetric  solution. 

For  smaller  values  of  the  interfacial  friction  coefficient,  the  increased  nonlinear¬ 
ity  led  to  low  frequency  variability  associated  with  changes  in  the  zonal  extent  of  the 
recirculating  gyres.  Formation  of  strong  rings  in  the  area  east  of  the  recirculations  and 
their  subsequent  westward  propagation  through  the  gyres  provided  intense  potential 
vorticity  mixing  across  the  front  and  in  the  area  of  the  recirculations.  Rings  were 
strong  enough  not  to  be  completely  dissipated  in  the  recirculations;  upon  reaching 
a  western  boundary,  they  rejoined  the  eastward  jet.  Therefore,  the  model  exhibited 
periods  of  large  recirculations  and  weak  eddy  activity,  when  the  model  trajectory  was 
in  the  vicinity  of  the  antisymmetric  fixed  point,  and  periods  of  small  recirculations 
and  strong  meander  and  ring  activity. 

With  further  decrease  in  the  interfacial  friction,  the  model  displayed  two  dis¬ 
tinct  energy  states.  A  high  energy  state  had  strong  variability  associated  with  well- 
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developed  recirculations  of  various  zonal  extents  described  in  the  above  paragraph. 
The  low  energy  state  corresponded  to  a  flow  regime  where  the  separation  point  of  the 
eastward  jet  exhibited  large  meridional  excursions.  Weak  recirculations  existed  west 
of  1000  km  only  in  the  time  mean  sense.  The  strong  eddy  field  was  dominated  by 
large  amplitude,  westward  propagating  meanders  whose  interaction  with  the  western 
boundary  current  produced  strong  meridionally  propagating  rings.  In  comparison, 
McCalpin  and  Haidvogel  [25]  identified  three  preferred  dynamical  regimes  each  with 
a  distinct  energy  level.  Although  the  high  and  the  medium  regime  resembled  the  vari¬ 
ability  associated  with  the  high  energy  state  of  the  regional  model,  the  low  energy 
state  was  quite  different  from  the  one  of  the  regional  model.  This  was  partly  because 
the  outflow  conditions  of  the  colliding  jets  model  forced  the  eastward  jet  to  stay 
coherent,  and  therefore  the  model  never  reached  the  point  when  the  jet  disappears 
almost  right  after  separation,  as  in  a  low  energy  state  of  the  wind-driven  model. 

For  low  values  of  the  interfacial  friction  coefficient,  the  flow  stayed  in  the 
dynamical  regime  described  as  a  low  energy  state  in  the  above  paragraph.  The 
amplitude  of  the  eddy  field  was  very  strong,  and  small  and  weak  recirculations  existed 
only  in  the  time  mean  sense. 

In  a  regional  colliding  jets  model,  the  role  of  wind  forcing  was  reduced  to 
prescribing  a  particular  distribution  of  streamfunction  and  potential  vorticity  on  the 
open  boundaries.  The  local  nature  of  such  forcing  allowed  isolating  the  role  of  internal 
barotropic  dynamics  of  the  jet  and  its  interaction  with  the  western  boundary  current 
as  mechanisms  responsible  for  the  existence  of  low  frequency  variability  and  multiple 
dynamical  regimes.  The  numerical  experiments  with  the  regional  model  showed  that 
low  frequency  variability  was  an  inherent  part  of  the  dynamics  where  nonlinearities 
were  not  too  weak  or  too  strong.  The  presence  of  the  western  wall  was  essential, 
since  no  low  frequency  variability  or  multiple  dynamical  regimes  were  observed  in  a 
barotropic,  zonally  asymmetric  model  of  Jayne  et  al  [18],  which  was  forced  by  an 
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unstable  zonal  jet  on  the  open  western  boundary  and  a  stable  zonal  jet  on  the  open 
eastern  boundary. 

One  of  the  most  important  properties  of  the  regional  model  was  a  constant 
energy  input,  which  was  given  by  the  boundary  integrated  energy  flux  vector.  In 
comparison,  in  the  wind-driven  models,  the  forcing  had  global  character,  and  the 
energy  input  strongly  depended  on  the  solution  form.  Scott  and  Straub  [34]  discussed 
a  correlation  mechanism  between  symmetric  wind-forcing  and  a  geostrophic  current 
which  led  to  ability  of  steady  antisymmetric  solutions  to  extract  more  energy  than 
nonsymmetric  solutions.  The  experiments  with  various  wind  asymmetry  parameter  in 
the  wind-driven  model  of  McCalpin  and  Haidvogel  [25]  showed  that  for  highly  asym¬ 
metric  wind  forcing,  the  model  was  not  able  to  reach  high  energy  states.  Since  the 
regional  model  did  not  have  any  mechanism  that  would  help  to  restore  the  symmetry 
of  the  solution,  it  was  harder  for  the  model  trajectory  to  reach  and  to  stay  in  the 
vicinity  of  the  antisymmetric  steady  solution. 

The  reduced-gravity  model  had  only  the  mechanism  of  barotropic  instability, 
and  therefore  the  question  remains  whether  the  presence  of  baroclinic  instability  could 
significantly  affect  low  frequency  variability.  The  next  Chapter  addresses  this  question 
in  a  framework  of  a  two  layer,  colliding  jets  model. 

In  the  present  study,  we  used  a  free  slip  boundary  condition  on  the  western 
boundary.  As  Berio ff  and  McWilliams  [3]  showed,  the  presence  of  a  no-slip  boundary 
condition  tends  to  destabilize  western  boundary  currents.  Eddies  produced  by  such 
currents  can  be  advected  towards  the  middle  of  the  domain  and  may  interact  with 
rings  that  propagate  towards  the  western  wall  from  the  model’s  interior.  Hence,  the 
nature  of  low  frequency  variability  and  of  multiple  dynamical  regimes  might  change 
with  the  use  of  no-slip  boundary  condition  due  to  more  complicated  dynamics  of  the 
western  boundary  current.  This  problem,  however,  was  not  addressed  in  the  present 
thesis. 
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Chapter  5 


Results  for  a  Two  Layer  Colliding 
Jets  Model 

5.1  Introduction 

In  this  Chapter,  we  present  the  results  of  the  numerical  experiments  with  a  two-layer, 
colliding  jets  model,  which  includes  the  mechanisms  of  both  barotropic  and  baroclinic 
instabilities.  The  question  is  how  the  latter  affects  the  existence  and  the  character  of 
low  frequency  variability  and  whether  multiple  dynamical  regimes  are  still  possible. 

As  we  showed  in  Chapter  2,  baroclinic  instability  is  more  efficient  than  barotropic 
instability  in  inducing  westward  flows  during  the  spin  down  of  a  quasi-geostrophic, 
two-layer,  zonally  symmetric  jet.  In  strongly  baroclinic  flows,  it  was  eddy  heat  fluxes 
that  drove  the  lower  layer  recirculations  through  the  residual  meridional  circulation. 
For  strongly  barotropic  flows,  the  mechanisms  were  nonlinear  eddy-eddy  and  eddy- 
mean  flow  interactions  that  produced  weak  westward  flows.  As  we  showed  in  the 
previous  Chapter,  the  dynamics  of  rings  was  also  crucial  for  the  existence  of  multi¬ 
ple  dynamical  regimes  in  the  reduced- gravity,  colliding  jets  model.  The  mechanism 
of  baroclinic  instability  associated  with  the  energy  transfer  from  the  vertical  shear 
to  the  eddy  field  might  be  considered  as  an  additional  “degree  of  freedom”  for  the 
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two-layer  colliding  jets  model  and  therefore  may  substantially  affect  the  character  of 
the  low  frequency  variability  and  the  existence  of  multiple  dynamical  regimes. 

The  formulation  of  the  model  and  its  governing  equations  are  presented  in  sec¬ 
tion  3.2.1.  The  dissipation  terms  are  given  by  V1  =  and  V2  =  -A2V2ip2  - 

AeVeip2,  where  As  is  a  biharmonic  diffusion  coefficient  and  A2  is  a  bottom  friction 
coefficient. 


5.2  Comments  on  Strong  Bottom  Friction 

As  we  shall  see,  bottom  friction  in  the  two-layer  model  physically  plays  an  opposite 
role  to  that  of  interfacial  friction  in  the  reduced-gravity  model;  as  bottom  friction 
increases,  the  upper  layer  flow  behaves  as  that  of  the  reduced-gravity  model  with  weak 
interfacial  friction.  The  physics  can  be  illustrated  by  the  following  scaling  argument. 
Suppose  that  bottom  friction  is  strong  in  the  two-layer  model,  and  therefore  the 
amplitude  of  the  lower  layer  streamfunction,  ^>2,  is  small.  Assume  also  that  the  lower 
layer  depth  is  large,  so  that  the  lower  layer  Froude  number  is  small.  With  the  above 
assumptions,  the  dominant  balance  in  the  lower  layer  potential  vorticity  equation  is 
between  the  time  changes  in  the  upper  layer  streamfunction  and  the  bottom  friction 
operator  applied  to  the  lower  layer  streamfunction: 

O  I 

F2~ - i2V^2-  (5.1) 

Rewriting  the  upper  layer  potential  vorticity  field  as  ITi  ~  /3y  +  V2V>i  -  Fifa,  the 
governing  equation  for  the  upper  layer  flow  becomes: 

5ITi 

-to+Jfa’  no  ~  -Avfyi  +  s1  +  ^-1{f1f2v-2J^}</,1.  (5.2) 

The  last  term  in  the  above  equation  appears  as  frictional  forcing  exerted  by  the  lower 
layer  flow  on  the  upper  layer  flow.  The  bottom  friction  coefficient  enters  with  the 
power  of  -1,  which  indicates  that  with  an  increase  in  A2,  the  upper  layer  flow  dy¬ 
namics  becomes  similar  to  that  of  a  reduced-gravity  model  with  nonzero  interfacial 
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friction.  This  explains  why  for  large  values  of  the  bottom  friction  coefficient,  the 
dynamical  regimes  of  a  time-dependent  two-layer  model,  described  in  the  following 
section,  closely  resemble  those  of  a  reduced-gravity  model  with  weak  interfacial  fric¬ 
tion. 


5.3  Sensitivity  Study  vs  Bottom  Friction  Coeffi¬ 
cient 

Case:  a  —  0 

The  behavior  of  the  two-layer  model  was  investigated  with  a  set  of  various  bottom 
friction  parameters.  Table  5.1  classifies  the  resulting  state  vs  A2  for  the  case  when  the 
lower  layer  inflow  velocity  was  chosen  zero,  so  that  the  model  was  forced  only  through 
the  upper  layer  inflow  and  outflow  conditions.  In  all  the  cases,  with  an  increase  in  the 
values  of  the  bottom  friction  coefficient,  the  model  showed  the  following  progression 
of  states: 

•  State  A.  The  model  has  a  unique  dynamical  regime  when  the  flow  consists  of 
the  large  upper  layer  recirculations  and  smaller,  eddy-driven  lower  layer  recir¬ 
culations.  The  total  energy  varies  by  only  3%,  which  indicates  weak  variability 
associated  with  changes  in  the  meridional  position  of  the  front  and  in  the  lon¬ 
gitudinal  extent  of  both  the  upper  and  the  lower  layer  recirculating  gyres.  The 
example  of  the  mean  streamfunction  field  is  given  in  Figure  5-5a,c.  The  detailed 
discussion  of  the  state  is  given  in  section  5.4. 

•  State  B.  The  upper  layer  flow  structure  consists  of  an  eastward  jet  and  two 
recirculating  gyres.  The  flow  in  the  lower  layer  is  weak  because  of  stronger  bot¬ 
tom  friction.  The  longitudinal  size  of  the  recirculations  changes  with  time  thus 
giving  rise  to  strong  low  frequency  variability.  The  latter  displays  itself  in  up  to 
16%  deviations  of  the  total  energy  from  its  mean  value.  The  amplitude  and  the 
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spatial  structure  of  the  eddy  field  are  strongly  time-dependent.  The  upper  layer 
eddy  field  is  characterized  by  the  eastward  and  westward  propagating  meanders. 
Their  interactions  give  birth  to  rings  in  the  area  east  of  the  recirculation  gyres 
in  the  upper  layer.  The  example  of  a  typical  time  mean  streamfunction  field  is 
given  in  Figure  5-1,  when  *42  =  5  •  10-6  sec-1. 

•  State  C.  The  flow  structure  displays  well-developed  low  frequency  variability 
associated  with  two  dynamical  regimes:  high  energy  state  B  or  low  energy  state 
D. 

•  State  D.  The  upper  layer  flow  is  characterized  by  an  eastward  jet  with  large 
meridional  excursions  in  the  separation  point.  The  corresponding  eddy  field 
is  dominated  by  large-amplitude,  westward  propagating  meanders  and  strong 
rings  forming  only  on  the  western  boundary.  The  recirculations  are  absent 
in  instantaneous  fields.  The  lower  layer  flow  is  characterized  by  two  counter¬ 
recirculations  that  are  very  weak  and  are  strongly  frictional.  The  example  of  a 
typical  time  mean  streamfunction  field  is  given  in  Figure  5-1,  when  ^42  =  5  TO-5 
sec-1. 

The  next  section  presents  the  detailed  analysis  of  a  reference  run  with  State 
A,  when  A2  =  10-8  sec-1.  In  states  B,  C  and  D,  the  lower  layer  flow  is  strongly 
frictional,  and  hence  these  cases  are  of  not  much  interest.  In  terms  of  the  upper  layer 
dynamics,  States  B ,  C  and  D  of  the  two-layer  model  correspond  to  States  D,  E,  F  of 
the  reduced- gravity  model.  Figure  5-3  shows  instantaneous  streamfunction  fields  for 
state  C,  when  A2  =  7.5  •  10-6  sec-1.  Clearly,  the  upper  layer  flow  structure  suffers 
significant  changes  with  time.  When  the  energy  is  high,  the  flow  is  characterized  by 
large  upper  layer  recirculations.  As  the  recirculations  decrease  in  size,  the  strength 
of  the  eddy  field  increases,  giving  rise  to  larger  amplitude  meanders  and  rings.  When 
the  energy  is  low,  the  flow  is  dominated  by  large  amplitude,  westward  propagating 
meanders  and  strong  rings  forming  on  the  western  boundary. 
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A  dim 

(sec  L) 

STATE 

1.0 

h—1 

O 

1 

00 

A 

2.0 

10"8 

A 

3.0 

10-8 

A 

5.0 

10-8 

A 

1.0 

10-7 

A 

8.0 

10-7 

B 

1.0 

10~6 

B 

3.0 

10“6 

B 

5.0 

10-6 

B 

7.5 

10~6 

C 

1.0 

10~5 

C 

3.0 

10"5 

D 

5.0 

10-5 

D 

Table  5.1:  Flow  character  for  various  dimensional  bottom  friction  coefficients,  when 
the  lower  layer  inflow  is  zero,  i.e.  a  =  0. 


Case:  a  —  0.1 

Table  5.2  classifies  the  flow  regimes  for  the  case  when  the  lower  layer  inflow  stream- 
function  was  10  times  weaker  than  the  corresponding  upper  layer  inflow  streamfunc- 
tion,  i.e.  a  =  0.1.  Figure  5-4  shows  time  mean  streamfunction  fields  for  both  layers. 
For  small  values  of  bottom  friction  parameter,  the  flow  was  characterized  by  large  re¬ 
circulations  in  the  upper  layer  and  small  recirculations  in  the  lower  layer  either  pressed 
against  the  western  boundary  currents  or  located  north  and  south  of  the  eastward 
jet.  There  was  not  much  low  frequency  variability  in  this  case.  For  larger  values  of 
the  bottom  friction  coefficient,  the  western  boundary  current  separated  earlier  from 
the  wall,  giving  rise  to  a  wide,  frictional  eastward  flowing  jet.  The  upper  layer  flow 
with  large  recirculations  showed  strong  low  frequency  variability  associated  with  two 
dynamical  regimes  as  described  earlier. 
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Afm  (sec”1) 

STATE 

1.0-  1CT8 

A 

1.0- 10~7 

A 

5.0  •  10-7 

C 

1.0-  io-6 

C 

5.0  -10~6 

c 

Table  5.2:  Flow  character  for  various  dimensional  bottom  friction  coefficients,  when 
the  amplitude  of  the  lower  layer  inflow  streamfunction  is  10%  of  that  of  the  upper 
layer  inflow  streamfunction,  i.e.  a  =  0.1. 


Steady  Solutions 

An  attempt  to  find  the  antisymmetric  steady  solutions  of  the  two-layer  model  by  using 
a  numerical  code  with  enforced  antisymmetry  failed.  The  resulting  statistically  steady 
solutions  were  characterized  by  recirculating  gyres  that  reached  the  sponge  layer  near 
the  eastern  wall.  The  total  energy  never  attained  a  steady  state  due  to  the  presence 
of  a  nonzero  antisymmetric  eddy  field.  In  the  original  time-dependent  model  without 
artificial  antisymmetry,  the  zonal  extent  of  the  recirculating  gyres  varied  between  500 
km  and  2500  km  depending  on  the  value  of  the  bottom  friction  coefficient. 

The  numerical  experiments  with  the  reduced-gravity,  colliding  jets  model  in 
Chapter  4  showed  that  biharmonic  diffusion  was  not  effective  in  dissipating  energy, 
and  therefore  interfacial  drag  was  required  to  reach  a  statistically  steady  solution.  In 
the  present  two-layer  model,  the  upper  layer  has  only  biharmonic  diffusion,  which  may 
explain  the  absence  of  a  steady  antisymmetric  solution  with  a  reasonable  zonal  extent 
of  the  recirculating  gyres.  Therefore,  the  presence  of  bottom  friction  in  the  lower  layer 
and  the  time-dependence  of  the  upper  layer  potential  vorticity  field  are  required  for 
the  existence  and  maintenance  of  the  observed  statistically  steady  state.  Perhaps,  the 
model  could  have  had  steady  antisymmetric  solutions  if  a  different  frictional  operator, 
such  as  lateral  diffusion,  were  applied. 
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It  is  unclear  whether  the  model  has  other,  asymmetric,  steady  states;  this 
problem  is  beyond  the  scope  of  the  present  thesis.  One  may  argue  that,  if  they 
exist,  the  model  might  display  multiple  dynamical  regimes.  However,  when  both 
barotropic  and  baroclinic  instabilities  are  present,  the  model  dynamics  becomes  more 
complicated,  and  the  steady  states  might  not  be  useful  in  describing  the  behavior 
of  system.  State  A  can  be  considered  as  a  strange  attractor  of  a  quasi-geostrophic, 
turbulent,  two-layer  model  with  weak  bottom  friction.  As  we  show  in  the  next  section, 
in  this  parameter  regime,  the  model  does  not  have  multiple  dynamical  states,  but 
instead  it  is  characterized  by  only  one  dynamical  regime  with  weak  low  frequency 
variability. 

5.4  Reference  Run 

5.4.1  Dimensional  Parameters 

The  reference  run  for  the  two  layer  model  has  the  following  dimensional  parameters: 

L  =  2400  km, 

Lx  =  4667  km, 

Ldef  =  47.64  km, 

A$m  =  8  •  1010  m4sec-1, 

Afm  =  1  •  10-8  sec-1, 

D\  =  600  m, 

D2  =  4000  m, 

U  =  1  m  sec-1. 

The  initialization  procedure  discussed  in  section  3.2.9  was  used  to  prescribe  the 
boundary  distributions  (3.4)  and  (3.5)  of  the  streamfunction  field  and  its  derivatives. 
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Figure  5-2:  Upper  layer  (a)  and  lower  layer  (b)  time  mean  streamfunction  field,  when 
a  =  0  and  A2  —  5  •  10-5  sec-1. 
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Figure  5-3:  Instantaneous  upper  layer  streamfunction  (left  column)  and  corresponding 
lower  layer  streamfunction  (right  column),  when  a  =  0  and  Ai  =  5  •  10~5  sec-1.  The 
structure  of  the  flow  changes  between  a  state  with  well  developed  recirculations  in  the 
upper  layer  and  a  weakly-meandering  jet  to  a  state  with  no  upper  layer  recirculations 
and  large  amplitude,  westward  propagating  meanders. 
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Figure  5-4:  Upper  layer  (a,c,e,g)  and  lower  layer  (b,d,f,h)  time  mean  streamfunction 
field,  when  A2  =  10-8  sec-1,  A2  =  10-7  sec-1,  A2  =  5  ■  10-7  sec-1,  and  A2  =  5  •  10-6 
sec-1  respectively.  Lower  layer  inflow  in  nonzero:  a  =  0.1. 


The  nondimensional  amplitude  of  the  inflow  was  chosen  as  Ri  =  1  for  the  upper  layer 
and  R2  =  0  for  the  lower  layer.  Hence,  the  energy  input  into  the  model  domain  was 
only  due  to  the  upper  layer  inflow,  and  the  main  dissipation  mechanism  was  provided 
by  bottom  friction.  Since  the  interfacial  drag  was  zero,  the  circulation  in  the  lower 
layer  was  induced  due  to  horizontal  transfer  of  heat  by  the  eddy  fluxes. 

The  model  was  spun  up  for  30  years,  until  a  statistically  steady  state  was 
reached.  An  additional  424  years  of  integration  allowed  collection  of  the  necessary 
statistics. 

5.4.2  Time  Mean  Flow 

As  shown  later,  the  model  has  one  preferred  energy  state  associated  with  a  unique 
dynamical  regime.  Figure  5-5  shows  the  time  mean  fields  of  streamfunction  and  po¬ 
tential  vorticity.  The  upper  layer  flow  is  characterized  by  two  symmetric  recirculating 
gyres,  which  are  2000  km  long  and  400  km  wide.  The  time  mean  potential  vorticity 
field  reveals  two  regions  of  anomalous  potential  vorticity  indicating  waters  brought 
by  the  western  boundary  currents  from  the  north  and  south  into  the  recirculation 
area.  The  upper  layer  recirculations  are  also  well-defined  in  the  instantaneous  fields. 

The  lower  layer  flow  has  a  more  complicated  spatial  structure.  There  are  two 
tight  recirculations,  2000  km  long  and  200  km  wide,  in  the  middle  of  the  domain. 
Two  weak  counter-rotating  gyres  are  located  to  the  north  and  south  of  the  recircu¬ 
lations.  The  time-mean  flows  in  the  lower  layer  are  driven  by  the  eddy  heat  fluxes, 
and  therefore  they  exist  only  in  the  time  mean  sense;  the  instantaneous  fields  are 
dominated  by  strong  eddies.  Outside  the  recirculation  regions,  the  lower  layer  poten¬ 
tial  vorticity  field  is  dominated  by  its  planetary  component,  whereas  in  the  area  of 
the  recirculations  it  is  homogenized.  Therefore,  near  the  edges  of  the  eastward  jet, 
the  time  mean  recirculations  have  the  same  sense  of  rotation  in  the  upper  and  lower 
layers,  whereas  further  south  and  north  the  time  mean  lower  layer  gyres  rotate  in 
the  opposite  direction  to  the  upper  layer  flows,  creating  a  region  of  strong  baroclinic 
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Figure  5-5:  Contours  of  upper  (a,b)  and  lower  (c,d)  layer  streamfunction  and  potential 
vorticity  fields.  Parameters  are  defined  in  section  5.4.1. 
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flow.  The  average  total  transport  by  both  recirculations  is  about  40  Sv,  partitioned 
between  16  Sv  in  the  lower  layer  and  24  Sv  in  the  upper  layer. 

5.4.3  Energetics 

The  histograms  (Figure  5-6)  based  on  the  time  series  of  kinetic,  potential  and  total 
energy  indicate  that  the  model  has  a  unique  energy  state.  To  see  whether  the  dis¬ 
tributions  are  normal,  a  MATLAB  QQPLOT  routine  was  used  to  plot  the  quantities 
of  a  given  time  series  versus  the  quantiles  of  the  generated  time  series  with  normal 
distribution  (Figure  5-7).  The  upper  layer  kinetic  energy  and  potential  energy  time 
series  are  close  to  normal  distribution.  The  tails  of  the  lower  layer  kinetic  energy 
time  series  do  not  fall  within  the  error  bars  and  therefore  indicate  significant  devia¬ 
tion  from  normal  distribution.  The  latter  may  be  the  result  of  the  presence  of  strong 
vortices  in  the  flow  structure. 

On  average,  total  energy  is  partitioned  between  87%  potential,  9%  upper  layer 
kinetic  and  4%  lower  layer  kinetic  energy.  Therefore,  the  time  series  of  total  and 
potential  energy  are  well  correlated.  Computed  by  using  a  MATLAB  SPECTRUM 
routine  with  a  hanning  window  of  20  and  100  years,  the  potential  energy  power  density 
spectrum  (Figure  5-8)  is  red  with  a  distinct  change  in  the  slope  at  approximately 
2  months.  There  is  a  distinct  peak  at  approximately  1.4  years.  Another  peak  is 
associated  with  a  time  scales  larger  than  10  years.  Hence,  similar  to  the  reduced- 
gravity  model,  the  two-layer  model  exhibits  low  frequency  variability.  However,  there 
is  only  one  dynamical  regime  in  the  model,  and  it  is  associated  with  large  recirculating 
gyres  in  the  upper  layer  and  with  smaller  eddy-driven  gyres  in  the  lower  layer. 

Peaks  in  the  upper  and  lower  layer  kinetic  energy  spectrum  (Figure  5-9)  are 
located  on  the  time  scale  of  1.4  years.  The  lower  layer  kinetic  energy  has  a  con¬ 
stant  slope  at  time  scales  less  than  1.4  years,  whereas  the  upper  layer  kinetic  energy 
spectrum  changes  the  slope  at  2  months.  The  latter  must  be  related  to  the  internal 
instability  time  scale. 
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Figure  5-6:  Histograms  of  upper  layer  kinetic  (a),  lower  layer  kinetic  (b),  potential 
(c)  and  total  energy  (d).  Parameters  are  defined  in  section  5.4.1. 
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Figure  5-7:  The  quantile-quantile  plot  for  various  time  series  (solid  lines)  vs  normal 
distribution,  (a):  upper  layer  kinetic  energy,  (b):  lower  layer  kinetic  energy,  (c): 
potential  energy.  Dashed  line  denotes  error  bars.  Parameters  are  defined  in  section 
5.4.1. 


Figure  5-10  summarizes  the  domain-integrated  time  mean  energy  conversion 
in  the  model.  The  energy  input  into  the  model  domain  is  due  to  the  upper  layer 
inflows.  The  dissipation  is  provided  through  biharmonic  diffusion  and  sponge  ab¬ 
sorption  in  both  layers  and  through  bottom  friction.  57%  of  energy  is  dissipated 
through  perturbations  and  42%  is  dissipated  by  the  mean  flow.  Total  perturbation 
energy  is  almost  equally  partitioned  between  ,  K2  and  P  .  The  mean  flow  energy 
is  partitioned  between  91%  of  potential  energy  and  8%  and  1%  of  upper  and  lower 
layer  kinetic  energy.  In  the  domain-integrated  sense,  the  barotropic  conversion  term 
is  two  orders  of  magnitude  larger  than  the  baroclinic  conversion  term.  The  bulk  of 
barotropic  conversion  occurs  in  the  upper  layer,  since  the  integral  of  BTC  in  the  lower 
layer  is  less  than  1%  of  the  total  value.  Bcrloff  and  McWilliams  [4],  in  their  two-layer 
wind-driven  model  with  symmetric  forcing  and  lateral  diffusion,  also  observed  strong 
barotropic  dominance,  when  the  flow  stayed  in  the  regime  of  strong  jet  penetration 
and  large  antisymmetric  recirculations.  In  the  wind-driven,  double-gyre  model  of 
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Figure  5-8:  Potential  energy  power  density  spectrum  (a)  and  frequency  times  poten¬ 
tial  power  density  spectrum  (c)  for  window  of  20  years.  Dashed  line  denotes  confi¬ 
dence  interval  at  95%.  (b)  and  (d):  same  for  a  100-year  long  window.  Parameters 
are  defined  in  section  5.4.1. 
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Figure  5-9:  Upper  layer  kinetic  energy  power  density  spectrum  (a)  and  frequency 
times  upper  layer  kinetic  power  density  spectrum  (c)  for  window  of  30  years.  Dashed 
line  denotes  confidence  interval  at  95%.  (b)  and  (d):  same  the  lower  layer.  Parameters 
are  defined  in  section  5.4.1. 


Figure  5-10:  Energy  conversion  diagram.  Parameters  are  defined  in  section  5.4.1 

Holland  [15],  the  baroclinic  conversion  term  was  also  two  orders  of  magnitude  smaller 
than  the  barotropic  conversion  one  in  the  experiment  with  small  lateral  diffusion. 

5.4.4  Time  Mean  Eddy  Fields 

Figure  5-11  presents  the  energy  density  of  the  mean  and  eddy  fields.  The  maxima  in 
the  upper  layer  mean  kinetic  energy  are  aligned  with  the  eastward  jet  west  of  2500 
km  and  with  the  colliding  western  boundary  currents.  The  lower  layer  mean  kinetic 
energy  has  three  maxima:  one  on  the  axis  of  the  jet  east  of  2000  km  and  the  other  two 
in  the  region  of  the  recirculating  gyres.  Two  pools  of  maximum  mean  potential  energy 
are  associated  with  the  recirculating  gyres.  Both  perturbation  potential  energy  and 
upper  layer  perturbation  kinetic  energy  have  maxima  on  the  axis  of  the  jet,  which 
reflects  the  importance  of  meanders  in  the  dynamics  of  the  flow.  The  lower  layer 
perturbation  kinetic  energy  has  two  maxima  in  the  western  part  of  the  recirculating 
gyres  and  along  the  western  wall. 
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Figure  5-12:  Energy  level  (in  %)  for  different  eigenmodes  resulted  from  an  extended 
EOF  analysis  for  a  streamfunction  fields,  where  time  scales  smaller  than  10  years 
were  filtered  out.  Dashed  line  denotes  a  95%  energy  level. 

5.4.5  EOF  Analysis  of  Low  Frequency  Variability 

To  understand  the  nature  of  low  frequency  variability  appearing  as  a  wide  peak  on 
the  time  scales  greater  than  10  years  in  the  potential  energy  spectra,  an  extended 
spatial  Empirical  Orthogonal  Function  (EOF)  decomposition  was  performed  for  a 
sub-domain  where  the  variability  was  maximum:  2917  km  x  1312  km  centered  near 
y  =  0.  The  streamfunction  fields  were  first  interpolated  on  to  a  regular  grid  of  32  x  32 
points  corresponding  to  the  above  sub-domain  (the  interpolation  grid  was  different 
from  the  model  grid).  Second,  a  finite  impulse  response  filter  was  applied  to  the 
time  series  to  exclude  time  scales  smaller  than  10  years.  Fourteen  spatial  EOFs  were 
required  to  reach  a  95%  energy  level  (Figure  5-12). 

The  first  EOF  is  symmetric  in  both  layers  with  a  variance  accounting  for  almost 
45%  of  the  total  energy  level.  The  corresponding  spatial  structure  of  the  mode  and 
its  amplitude  are  shown  in  Figures  5-13a-c.  The  mode  structure  indicates  variability 
associated  with  either  northward  or  southward  shifts  in  the  position  of  the  upper  layer 
jet  axis  in  the  area  between  the  recirculating  gyres.  The  position  of  the  separation 
point  is  fixed  in  the  middle  of  the  domain,  which  leads  to  the  jet  flowing  either  in  the 
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north-eastward  or  in  the  south-eastward  directions  after  separation  from  the  western 
boundary.  In  the  lower  layer,  the  mode  tends  to  perturb  the  front  near  the  separation 
point  and  to  strengthen  either  the  southern  or  northern  recirculations.  The  spectrum 
of  the  time-dependent  amplitude  indicates  a  wide  maximum  on  the  time  scales  larger 
than  10  years  with  a  slow  decay  towards  smaller  frequencies.  Longer  time  series  is 
required  to  isolate  the  variability  on  the  scale  close  to  100  years. 

The  other  50%  of  the  95%-energy  level  are  partitioned  between  thirteen  EOFs, 
with  the  second  and  the  third  one  accounting  for  10.6%  and  10.3%  respectively  (Fig¬ 
ures  5-13d-i).  The  spatial  structure  of  the  second  mode  indicates  variability  in  the 
area  where  the  jet  splits;  it  is  associated  with  either  southward  and  northward  excur¬ 
sion  of  the  front  and  with  “wrapping”  of  the  recirculations  around  each  other.  In  the 
lower  layer,  the  second  mode  tends  to  increase  the  longitudinal  extent  of  one  of  the 
recirculations  and  decrease  that  of  the  other  one.  Also,  it  perturbs  the  front  near  the 
separation  point. 

When  the  amplitude  of  the  third  EOF-mode  is  positive,  the  upper  layer  re¬ 
circulations  tend  to  increase  in  their  longitudinal  extent,  and  the  front  has  a  weak 
north-west  to  south-east  slope.  In  the  lower  layer,  the  symmetry  of  the  flow  breaks, 
since  western  part  of  the  northern  recirculation  and  the  eastern  part  of  the  southern 
one  tend  to  weaken.  When  the  amplitude  is  negative,  the  upper  layer  recirculations 
shrink  in  zonal  direction,  and  the  slope  of  the  front  is  from  south-west  to  north-east. 
The  lower  layer  recirculations  appears  weakened  in  the  eastern  part  for  the  northern 
one  and  in  the  western  part  for  the  southern  part. 

A  similar  EOF  decomposition  based  on  a  streamfunction  time-series  with  vari¬ 
ability  larger  than  20  years  filtered  out  require  ten  EOFs  for  the  95%  total  energy 
level,  with  the  most  energetic  mode  having  a  58%  energy  level.  The  spatial  structure 
and  the  amplitude  of  the  most  important  EOFs  compare  well  with  those  based  on  a 
10-year  low-pass  filter. 
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Figure  5-13:  The  time-dependent  amplitude  (a)  and  the  upper  (b)  and  the  lower 
layer  (c)  spatial  structures  of  the  first  EOF.  Same  for  the  second  (d-f)  and  third  (g-i) 
EOFs.  The  EOF  analysis  is  based  on  the  streamfunction  field  with  variability  larger 
than  10  years  filtered  out.  i  qq 


5.4.6  EOF  Analysis  of  High  Frequency  Variability 

Inter-annual  Variability 

To  study  the  inter-annual  variability  of  the  model,  we  performed  an  EOF  analysis  for 
the  following  two  streamfunction  time  series.  One  was  band-passed  filtered  for  the 
variability  between  1  and  10  years  and  the  other  one  for  the  variability  between  1 
and  2  years.  In  both  cases,  the  most  important  EOFs  were  remarkably  similar,  and 
the  number  of  modes  required  for  the  95%  energy  level  was  very  close:  28  and  25 
respectively.  Therefore,  most  of  the  inter-annual  variability  is  on  the  scales  between 
1  and  2  years,  which  is  supported  by  the  energy  spectra  plots  of  section  5.4.3.  Thus, 
here  we  present  the  results  of  only  the  second  EOF  decomposition. 

Figure  5-14  shows  variance  associated  with  25  EOFs  required  for  the  95% 
energy  level  with  the  first  three  functions  having  17%,  14%  and  11%  energy  levels 
respectively  (for  comparison:  for  the  time  series  filtered  between  1  and  10  years,  the 
first  three  EOFs  had  11%,  10%  and  8%  energy  levels).  The  spectra  of  the  time- 
dependent  amplitude  and  the  corresponding  spatial  structure  are  presented  in  Figure 
5-15.  All  modes  are  clearly  associated  with  the  meanders  of  approximately  1000  km 
length  scale.  The  peak  in  amplitude  spectra  is  around  1.3  years. 

Intra-annual  Variability 

To  study  the  inter-annual  variability  of  the  model,  we  performed  an  EOF  analysis  for 
the  following  two  streamfunction  time  series.  One  was  band-passed  filtered  for  the 
variability  between  1  and  12  months  and  the  other  one  for  the  variability  between  5 
and  12  months.  In  both  cases,  the  most  important  EOFs  were  remarkably  similar. 
Here  we  present  the  results  of  only  the  second  EOF  decomposition. 

Total  of  35  modes  required  to  reach  the  95%  energy  level.  The  first  three 
EOFs  accounted  for  24%,  18%  and  9%  of  the  total  energy  level.  Their  amplitude  and 
spatial  structures  are  presented  in  Figure  5-17.  The  modes  are  symmetric  wave-like 
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Figure  5-14:  Energy  level  (in  %)  for  different  eigenmodes  resulted  from  an  extended 
EOF  analysis  for  a  streamfunction  with  variability  band-passed  between  1  and  2 
years.  Dashed  line  denotes  a  95%  energy  level. 

structures  sitting  on  the  axis  of  the  upper  and  the  lower  layer  jet.  The  length  scale 
of  the  modes,  decaying  west  of  2000  km,  is  approximately  1000  km.  The  frequency 
multiplied  by  the  power  density  spectrum  of  the  time-dependent  amplitude  shows  a 
peak  at  8-9  months  for  all  three  modes.  Therefore,  the  first  three  modes  account 
for  51%  of  energy  and  are  responsible  for  the  variability  on  the  time  scales  around  9 
months.  The  other  modes  are  each  less  energetic,  but  together  they  account  for  the 
remaining  44%  of  the  total  energy  level. 

5.5  The  Role  of  Baroclinic  Instability  in  Low  Fre¬ 
quency  Variability 

5.5.1  Method 

To  study  the  role  of  baroclinic  instability  in  low  frequency  variability,  we  employed  a 
statistical  method  based  on  the  approach  of  Davis  and  Emanuel  [7].  Applied  to  the 
time  series  of  the  streamfunction  fields,  the  method  allows  to  define  zones  of  baroclinic 
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Figure  5-15:  The  time-dependent  amplitude  (a)  and  the  upper  (b)  and  the  lower  layer 
(c)  spatial  structures  of  the  first  EOF.  Same  for  the  second  (d-f)  and  third  (g-i)  EOFs. 
The  EOF  analysis  is  based  on  the  streamfunction  field  with  variability  between  1  and 
2  years.  202 
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Figure  5-16:  Energy  level  (in  %)  for  different  eigenmodes  resulted  from  an  extended 
EOF  analysis  for  a  streamfunction  with  variability  between  5  and  12  months.  Dashed 
line  denotes  a  95%  energy  level. 

growth  in  the  time  average  sense.  At  each  time  step,  the  streamfunction  and  potential 
vorticity  fields  are  decomposed  into  a  time  mean  flow  and  perturbations: 

ipi{x,y,t)  =  ^i(x,y)  +ipi(x,y,t)  (5.3) 

and 

Qi(x,y,t)  =  <ii{x,y)  +qi(x,y,t),  i  =  1,2.  (5.4) 

The  perturbation  streamfunction  is  further  rewritten  as  a  sum  of  a  field  associated 
with  the  upper  layer  potential  vorticity  anomaly  and  a  field  associated  with  the  lower 
layer  potential  vorticity  anomaly,  i.e. 

ip'i  =  +  fpi2\  *  =  !» 2 

where  is  the  solution  of  the  problem: 

'  (V2-  Wi1)+Fi#=9i, 

*  (V2-E2)#+F2^1)  =  0, 

and  V,|2'>  is  the  solution  of  the  problem: 

’  (V2  -  F\)vf }  =  0, 

k  (V2  -  F2)42)  +  =  q2. 


(5.5) 

(5.6) 

(5.7) 
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Figure  5-17:  The  frequency  multiplied  by  the  power  density  spectrum  of  the  square 
of  the  time-dependent  amplitude  (a)  and  the  contours  of  the  upper  (b)  and  the  lower 
layer  spatial  structures  (c)  of  the  first  EOF.  Same  for  the  second  (d-f)  and  third  (g-i) 
EOFs.  The  EOF  analysis  is  based  on  the  streamfunction  time  series  with  variability 
between  5  and  12  months.  The  spatial  structure  shown  for  the  sub-domain,  units  are 
in  km. 
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The  linearized  perturbation  enstrophy  equations  then  can  be  rewritten  as: 


where  j  is  a  geostrophic  velocity  field.  Zones  of  baroclinic 

growth  are  defined  by  ( x ,  y)  where  both  of  the  following  constraints  are  satisfied: 


(qW?  ■  V^j)  <  0, 
>  (?2u21}  •  V?2)  <  0, 


(5.10) 


where  (•)  denotes  time  average.  These  inequalities  can  be  interpreted  as  conditions  for 
an  increase  in  perturbation  enstrophy  due  to  a  down- gradient  eddy  potential  vorticity 
flux  produced  by  the  potential  vorticity  anomaly  in  the  other  layer. 


5.5.2  Applications 

The  above  method  was  applied  for  two  different  values  of  the  bottom  friction  coeffi¬ 
cient:  A2  =  10-8  sec-1  and  A2  —  3  •  10-6  sec-1.  Figure  5-18a  shows  that  in  the 
first  case,  there  are  two  regions  of  baroclinic  growth  located  in  the  area  of  westward 
recirculations  accompanied  by  other  two  small  zones  north  and  south  of  where  the 
jet  splits.  As  was  shown  in  the  previous  sections,  for  a  given  value  of  bottom  friction 
coefficient,  the  total  energy  varies  by  less  than  3%  from  its  mean  value,  and  the  low 
frequency  variability,  although  present,  is  weak. 

For  A2  =  3  •  10-6  sec-1,  small  zones  of  baroclinic  growth  are  located  on  the 
most  western  and  eastern  edges  of  the  recirculations,  where  the  flow  becomes  pre¬ 
dominantly  meridional.  The  total  energy  level  varies  by  16%,  and  the  low  frequency 
variability  is  strong. 

Therefore,  for  small  values  of  the  bottom  friction  parameter,  the  presence  of 
baroclinic  instability  in  the  recirculation  area  reduces  the  strength  of  low  frequency 
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variability.  In  these  experiments,  multiple  dynamical  regimes  do  not  occur  when 
baroclinic  instability  is  present  in  the  recirculation  region. 

5.6  Summary  and  Discussion 

In  Chapter  5,  we  described  the  results  of  the  numerical  experiments  with  a  two-layer, 
colliding  jets  model  which  showed  strong  dependence  of  the  flow  structure  on  the 
bottom  friction  parameter.  For  small  bottom  friction,  which  implied  a  physically  rel¬ 
evant  lower  layer  flow,  the  model  stayed  in  a  unique  dynamical  regime  associated  with 
the  well-developed  recirculations  in  both  layers.  To  understand  the  nature  of  low  and 
high  frequency  variability  observed  in  the  model,  an  extended  spatial  EOF  analysis 
was  performed  for  filtered  streamfunction  fields.  The  low  frequency  spatial  EOFs 
were  found  to  explain  decadal  and  longer  variability  in  the  position  of  the  eastward 
jet  axis,  as  well  as  weak  wrapping  and  zonal  pulsations  of  the  recirculation  gyres.  For 
quasi-annual  and  mesoscale  variability,  the  EOF  analysis  showed  the  importance  of 
wave-like  structures  sitting  on  the  axis  of  the  jet  between  the  recirculating  gyres  with 
the  length  scale  of  1000  km  and  the  time  scale  of  8  to  9  months. 

For  large  values  of  the  bottom  friction  coefficient,  the  lower  layer  flow  became 
strongly  frictional,  and  the  recirculations  disappeared.  However,  the  dynamics  of 
the  upper  layer  flow  became  more  complicated  displaying  both  multiple  dynamical 
regimes  and  low  frequency  variability  which  closely  resembled  those  of  the  reduced- 
gravity  model  with  weak  interfacial  friction.  Such  similarity  between  the  reduced- 
gravity  and  the  two-layer  models  was  possible  because  of  the  forcing  provided  by  the 
lower  layer  flow  on  the  upper  layer  flow  and  associated  with  strong  bottom  friction. 

The  results  of  the  two-layer  calculations  showed  that  for  physically  feasible 
values  of  the  bottom  friction  coefficients,  the  presence  of  baroclinic  instability  in 
the  recirculation  area  prohibited  the  existence  of  multiple  dynamical  regimes.  How¬ 
ever,  weak  low  frequency  variability  associated  with  changes  in  the  structure  and  the 
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(a) 


Figure  5-18:  Zones  of  baroclinic  growth  for  A2  —  10  8  see  1  (a)  and  A2  —  3-10  6 
sec-1  (b).  White  contours  show  upper  layer  time  mean  streamfunction  field. 
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strength  of  the  recirculating  gyres  was  still  possible.  In  contrast  with  the  reduced- 
gravity  case,  where  steady  antisymmetric  solutions  helped  to  understand  the  exis¬ 
tence  and  the  structure  of  a  high  energy  state  and  the  associated  dynamical  regime, 
an  attempt  to  find  steady  antisymmetric  solutions  by  employing  a  two-layer  model 
with  enforced  antisymmetry  failed.  In  disagreement  with  the  original  time  dependent 
model,  where  the  zonal  extent  of  the  recirculations  varied  between  500  km  and  2500 
km,  the  recirculations  reached  all  the  way  to  the  eastern  sponge  layer,  and  the  flow 
was  characterized  by  a  statistically  steady  nonzero  antisymmetric  eddy  field.  There¬ 
fore,  the  total  energy  level  never  reached  a  steady  state.  One  of  the  reasons  why 
the  antisymmetric  steady  state  did  not  exist  may  be  that  the  upper  layer  dissipation 
mechanism  was  presented  only  by  biharmonic  diffusion,  which  is  not  efficient  at  dis¬ 
sipating  energy.  Therefore,  the  balance  between  the  nonlinear  terms  and  dissipation 
was  not  possible,  and  time  variability  was  essential.  The  use  of  lateral  diffusion  might 
have  resulted  in  existence  of  steady  antisymmetric  solutions.  However,  the  question 
still  remains  whether  the  corresponding  steady  states  would  have  been  relevant  to  the 
behavior  of  the  time-dependent  model.  As  the  system  becomes  progressively  more 
complex,  as  in  the  case  with  both  barotropic  and  baroclinic  instability  at  work,  the 
model  might  become  purely  turbulent,  and  the  resulting  dynamical  state  may  be  a 
strange  attractor  for  a  given  parameter  regime. 

In  contrast,  in  the  two-layer,  wind-driven,  quasi-geostrophic  model  by  Berio ff 
and  McWilliams  [4],  rare  transitions  between  three  distinct  dynamical  states  were 
possible.  Compared  with  a  corresponding  reduced-gravity,  wind-driven  model,  the 
presence  of  baroclinic  instability  weakened  and  regularized  low  frequency  variability. 
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Chapter  6 


Conclusions 


This  thesis  studies  the  problems  of  generation  and  maintenance  of  recirculations  by 
the  Gulf  Stream  instabilities  using  two  approaches.  The  first  examines  the  develop¬ 
ment  of  recirculations  during  the  spin  down  of  a  zonally  symmetric,  unstable,  quasi- 
geostrophic  jet.  The  second  studies  the  problem  of  existence  and  maintenance  of  low 
frequency  variability  and  of  multiple  dynamical  regimes  due  to  the  mechanisms  of 
barotropic  and  baroclinic  instabilities. 

Chapter  2  addresses  the  first  problem  in  a  framework  of  the  reduced-gravity 
and  the  two-layer  doubly-periodic,  quasi-geostrophic  models.  The  recirculations  are 
robust  features  of  a  two-layer  model,  which  includes  the  mechanisms  of  both  baroclinic 
and  barotropic  instabilities.  The  relative  strength  of  the  recirculations,  defined  as 
the  ratio  of  the  total  westward  to  the  total  eastward  transport,  is  a  function  of  the 
nondimensional  parameters  of  the  model.  It  decreases  as  the  value  of  nondimensional 
p  increases.  This  decrease  in  the  recirculation  strengths  is  due  to  a  decrease  in 
the  zone  of  homogenized  potential  vorticity,  which  in  turn  is  due  to  (approximate) 
enstrophy  conservation. 

When  P,  7  and  total  energy  level  are  fixed,  the  strength  of  the  recirculations  is 
maximum  for  strongly  baroclinic  flows  and  minimum  for  intermediate  values  of  a  baro¬ 
clinic  parameter,  ■F  u12+f2Ui  >  where  iq  is  i— th  layer  nondimensional  velocity  amplitude 
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and  Fi  is  i-th  layer  Froude  number.  In  strongly  baroclinic  cases,  eddies  marginally 
equilibrate  the  flow  and  create  a  wide  zone  of  homogenized  potential  vorticity  in  the 
lower  layer.  This  result  is  in  agreement  with  the  study  of  a  weakly-nonlinear  equili¬ 
bration  of  Phillips’  model  with  weak  dissipation  (Pedlosky  [28]).  In  the  upper  layer,  a 
positive  potential  vorticity  gradient  on  the  axis  of  the  jet  decreases  by  less  than  10%, 
thus  appearing  as  an  efficient  barrier  to  mixing.  Most  vigorous  mixing  occurs  on  the 
edges  of  the  jet  in  the  upper  layer  and  under  the  axis  of  the  jet  in  the  lower  layer. 

As  the  initial  basic  jet  becomes  more  barotropic,  a  positive  potential  vorticity 
gradient  develops  on  the  axis  of  the  jet  in  the  lower  layer,  reflecting  a  more  stable 
flow  structure.  In  strongly  barotropic  jets,  there  is  vigorous  mixing  both  across  the 
positive  potential  vorticity  gradient  and  on  the  edges  of  the  jet,  leading  to  a  decrease 
by  more  than  45%  in  the  value  of  the  upper  layer  potential  vorticity  gradient  on  the 
axis  of  the  jet. 

The  linear  stability  analysis  applied  to  the  two-layer  jets  shows  that  in  strongly 
baroclinic  cases,  the  lower  layer  recirculations  are  driven  by  the  eddy  heat  fluxes 
through  the  residual  meridional  circulation.  In  the  upper  layer,  the  divergence  of  eddy 
potential  vorticity  fluxes  tends  to  simply  “diffuse”  the  jet  during  the  stage  of  linear 
growth.  Nonlinear  eddy-eddy  and  eddy-mean  flow  interactions  are  further  important 
in  inducing  the  upper  layer  recirculations.  The  linear  analysis  shows  that  the  strength 
of  the  residual  meridional  circulations  increases  with  an  increase  in  the  width  of  basic 
flow.  That  is  why  a  “mean  flow”  model  ( Flierl ,  personal  communication)  successfully 
reproduces  recirculations  when  applied  to  strongly  baroclinic  basic  flows. 

For  strongly  barotropic  two-layer  basic  flows,  which  are  strongly  supercritical, 
the  linear  eddy  potential  vorticity  fluxes  tend  to  diffuse  the  jet  in  both  layers  without 
any  tendency  for  the  formation  of  westward  flows  on  the  flanks.  As  the  basic  flow 
becomes  wider,  it  stabilizes.  Therefore,  strongly  nonlinear  regimes,  which  allow  for 
both  eddy-mean  flow  and  eddy-eddy  interactions,  are  required  for  the  emergence  of 
recirculations  for  strongly  barotropic  basic  flows. 
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Thus,  although  the  recirculations  develop  in  a  two-layer  model  for  both  barotropic 
and  baroclinic  flows,  the  processes  leading  to  their  generation  are  different.  In  baro- 
clinic  case,  eddy-mean  flow  interaction  is  powerful  enough  to  drive  the  recirculations. 

In  barotropic  cases,  strongly  nonlinear  eddy-mean  flow  and  eddy-eddy  interactions  are 
required  for  the  development  of  recirculations.  The  existence  of  two  different  mecha¬ 
nisms  is  responsible  for  a  non-monotonic  dependence  of  the  recirculation  strength  on 
the  baroclinic  parameter  in  the  nonlinear  experiments. 

In  the  reduced-gravity  model,  which  has  the  mechanism  of  barotropic  instabil¬ 
ity  only,  recirculations  form  as  a  result  of  a  nonlinear  spin  down  of  an  unstable  basic 
flow  only  for  the  case  of  an  infinite  dimensional  deformation  radius.  For  deformation 
radii  smaller  than  65  km,  the  total  westward  transport  is  less  than  5%  of  the  total 
eastward  transport. 

The  results  of  the  linear  stability  analysis  are  misleading  since  they  imply  a 
weak  tendency  for  the  formation  of  the  recirculations  for  small  dimensional  deforma¬ 
tion  radii  and  the  absence  of  the  recirculations  for  an  infinite  dimensional  deformation 
radius.  The  linear  stability  analysis  is  not  helpful  since  the  basic  flow  is  only  slightly 
supercritical  for  large  7,  implying  a  weak  eddy  field. 

The  second  part  of  the  thesis  examines  the  problem  of  existence  and  mainte¬ 
nance  of  low  frequency  variability  and  of  multiple  dynamical  regimes  as  the  result  of 
barotropic  and  baroclinic  instabilities. 

Chapter  3  presents  the  formulation  of  a  quasi-geostrophic,  colliding  jets  model 
which  is  forced  by  prescribed  inflow  and  outflow  through  the  open  boundaries.  The 
geometry  and  the  amplitude  of  such  boundary  forcing  is  important  for  the  resulting 
structure  of  a  time-dependent  solution.  The  inflow  consists  of  two  thin  jets  enter¬ 
ing  the  domain  from  the  north  and  the  south  and  flowing  along  the  closed  western 
boundary.  After  colliding  in  the  middle  of  the  domain,  the  western  boundary  cur¬ 
rents  give  rise  to  an  eastward  flowing  jet  that  leaves  the  domain  through  the  eastern 
boundary.  A  free  slip  boundary  condition  is  applied  on  the  western  boundary,  which 
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tends  to  suppress  local  instabilities  of  the  western  boundary  currents.  The  presence 
of  the  western  wall  and  western  boundary  currents  is  crucial  for  the  the  existence  of 
multiple  dynamical  regimes. 

The  model  has  two  types  of  friction:  biharmonic  diffusion  and  interfacial  drag 
in  a  reduced-gravity  formulation  and  biharmonic  diffusion  and  bottom  drag  in  a  two- 
layer  formulation.  In  addition,  the  sponge  layers  located  along  the  open  boundaries 
allow  the  absorption  of  the  outgoing  radiation  and  provide  a  smooth  transition  be¬ 
tween  the  interior  solution  and  the  boundary  conditions. 

Although  negligible  in  the  local  potential  vorticity  balance  of  the  recircula¬ 
tion  regions,  the  wind  is  nevertheless  important  for  the  existence  and  maintenance 
of  multiple  dynamical  regimes  in  the  wind-driven  models.  First,  the  model  energy 
input,  defined  as  the  domain-integrated  streamfunction  multiplied  by  the  wind-stress 
curl,  is  time-dependent.  When  the  nonlinear  interactions  are  at  work,  the  correlation 
between  the  spatial  structures  of  the  wind-stress  curl  and  the  solution  can  help  the 
model  trajectory  to  reach  a  linear  neighborhood  of  an  antisymmetric  fixed  point  and 
therefore  to  lead  the  model  towards  a  high  energy  state.  Second,  the  potential  vortic¬ 
ity  input  into  various  water  masses,  such  as  the  recirculations,  depends  strongly  on  the 
structure  of  the  solution.  For  nonsymmetric  solutions,  where  one  of  the  recirculations 
is  wrapping  around  the  other  one  and  is  crossing  the  line  of  zero  wind-stress  curl,  the 
potential  vorticity  forcing  integrated  over  a  closed  streamline  decreases  in  absolute 
values  compared  to  that  of  the  antisymmetric  solutions.  Thus,  in  the  wind-driven 
models,  there  are  two  possible  sources  of  low  frequency  variability:  wind  forcing  and 
instability  mechanisms. 

The  colliding  jets  model  is  highly  idealized,  however,  it  isolates  the  important 
mechanisms  of  instability:  barotropic  in  a  reduced-gravity  case  and  both  barotropic 
and  baroclinic  in  a  two-layer  case.  The  model  differs  from  the  wind-driven  ones  in  that 
the  energy  input  due  to  inflow  through  the  open  boundaries  is  time-independent.  The 
boundary  forcing  does  not  act  directly  in  the  recirculation  region  and  has  no  mecha- 
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nisms  that  would  locally  restore  the  symmetry  of  the  solution  structure.  Therefore, 
the  observed  high  energy  states  are  more  transient  than  those  in  the  wind-driven 
cases  with  symmetric  wind. 

Chapter  4  reports  the  results  of  the  numerical  experiments  with  the  reduced- 
gravity,  colliding  jets  model.  For  a  wide  range  of  interfacial  friction  coefficients, 
steady  antisymmetric  solutions  of  the  reduced-gravity,  colliding  jets  model  consist  of 
an  eastward  jet  flanked  by  two  recirculating  gyres.  For  a  set  of  chosen  parameters, 
the  zonal  extent  of  recirculations  and  the  total  energy  level  increase  with  a  decrease 
in  the  interfacial  friction  coefficient. 

The  initial  value  technique  was  used  to  calculate  the  exponentially  growing 
solutions  of  the  antisymmetric  states.  The  amplitude  of  the  unstable  mode  closely 
resembles  that  of  the  unstable  solution  of  the  wind-driven  reduced-gravity  model  of 
Primeau  [29]  and  consists  of  four  strong  cells  sitting  on  the  axis  of  the  jet  between 
the  recirculating  gyres  and  weaker  cells  in  the  recirculation  regions.  The  e-folding 
time  varies  between  4.8  years  and  417  days.  The  linear  growth  in  the  perturbation 
amplitude  halts  when  the  nonlinear  interactions  of  the  cells  become  strong.  These 
interactions  lead  to  formation  of  rings  in  the  area  where  the  jet  widens. 

Fully  nonlinear  calculations  show  that  with  an  increase  in  the  value  of  the 
interfacial  friction  coefficient,  the  model  becomes  weakly  unstable  and  displays  a 
unique  dynamical  regime  corresponding  to  an  eastward  jet  with  two  recirculating 
gyres.  For  near  critical  values  of  the  friction  coefficient,  the  eddy  field  is  dominated  by 
eastward  and  westward  propagating  meanders  of  small  amplitudes.  Rings  are  weak 
and  therefore  dissipate  during  their  propagation  through  the  recirculation  regions. 
The  energy  level  and  the  flow  structure  are  close  to  those  of  the  steady  antisymmetric 
solution. 

For  smaller  values  of  the  interfacial  friction  coefficient,  the  increased  nonlin¬ 
earity  leads  to  low  frequency  variability  associated  with  changes  in  the  zonal  extent 
of  the  recirculating  gyres.  Formation  of  strong  rings  in  the  area  east  of  the  recircula- 
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tions  and  their  subsequent  westward  propagation  through  the  gyres  provides  intense 
potential  vorticity  mixing  across  the  front  and  in  the  area  of  the  recirculations.  Rings 
are  strong  enough  not  to  be  completely  dissipated  in  the  recirculations;  upon  reach¬ 
ing  a  western  boundary,  they  rejoin  the  eastward  jet.  Therefore,  the  model  exhibits 
periods  of  large  recirculations  and  weak  eddy  activity,  when  the  model  trajectory  is 
in  the  vicinity  of  the  antisymmetric  fixed  point,  and  periods  of  small  recirculations 
and  strong  meander  and  ring  activity. 

With  further  decrease  in  the  interfacial  friction,  the  model  displays  two  distinct 
energy  states.  The  high  energy  state  has  strong  variability  associated  with  well- 
developed  recirculations  of  various  zonal  extents,  described  in  the  above  paragraph. 
The  low  energy  state  corresponds  to  a  flow  regime  where  the  separation  point  of 
the  eastward  jet  exhibits  large  meridional  excursions.  Weak  recirculations  exist  west 
of  1000  km  only  in  the  time  mean  sense.  The  strong  eddy  field  is  dominated  by 
large  amplitude,  westward  propagating  meanders  whose  interaction  with  the  western 
boundary  current  produces  strong  meridionally  propagating  rings.  In  comparison, 
McCalpin  and  Haidvogel  [25]  identified  three  preferred  dynamical  regimes  each  with 
a  distinct  energy  level.  Although  the  high  and  the  medium  regimes  resemble  the 
variability  associated  with  the  high  energy  state  of  the  regional  model,  the  low  energy 
state  is  quite  different  from  the  one  of  the  regional  model.  This  is  partly  because  the 
outflow  conditions  of  the  colliding  jets  model  force  the  eastward  jet  to  stay  coherent, 
and  therefore  the  model  never  reaches  the  point  when  the  jet  disappears  almost  right 
after  separation,  as  in  a  low  energy  state  of  the  wind-driven  model. 

For  low  values  of  the  interfacial  friction  coefficient,  the  flow  stays  in  the  dynam¬ 
ical  regime  described  as  a  low  energy  state  in  the  above  paragraph.  The  amplitude 
of  the  eddy  field  is  very  strong,  and  small  and  weak  recirculations  exist  only  in  the 
time  mean  sense. 

In  a  regional  colliding  jets  model,  the  role  of  wind  forcing  is  reduced  to  pre¬ 
scribing  a  particular  distribution  of  streamfunction  and  potential  vorticity  on  the 
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open  boundaries.  The  local  nature  of  such  forcing  allows  isolating  the  role  of  internal 
barotropic  dynamics  of  the  jet  and  its  interaction  with  the  western  boundary  current 
as  mechanisms  responsible  for  the  existence  of  low  frequency  variability  and  multiple 
dynamical  regimes.  The  numerical  experiments  with  the  regional  model  show  that 
low  frequency  variability  is  an  inherent  part  of  the  barotropic  dynamics  where  nonlin¬ 
earities  are  not  too  weak  or  too  strong.  The  presence  of  the  western  wall  is  essential, 
since  no  low  frequency  variability  or  multiple  dynamical  regimes  were  observed  in  a 
barotropic,  zonally  asymmetric  model  of  Jayne  et  al.  [18],  which  was  forced  by  an 
unstable  zonal  jet  on  the  open  western  boundary  and  a  stable  zonal  jet  on  the  open 
eastern  boundary. 

One  of  the  most  important  properties  of  the  regional  model  is  a  constant 
energy  input,  given  by  the  boundary  integrated  energy  flux  vector.  In  comparison, 
in  the  wind-driven  models,  the  forcing  has  global  character,  and  the  energy  input 
strongly  depends  on  the  solution  form.  Scott  and  Straub  [34]  discussed  a  correlation 
mechanism  between  symmetric  wind-forcing  and  a  geostrophic  current  which  led  to 
ability  of  steady  antisymmetric  solutions  to  extract  more  energy  than  nonsymmetric 
solutions.  The  experiments  with  various  wind  asymmetry  parameter  in  the  model 
of  McCalpin  and  Haidvogel  [25]  showed  that  for  highly  asymmetric  wind  forcing,  the 
model  was  not  able  to  reach  high  energy  states.  Since  the  regional  model  does  not 
have  any  mechanism  that  restores  the  symmetry  of  the  solution,  it  is  harder  for  the 
model  trajectory  to  reach  and  to  stay  in  the  vicinity  of  the  antisymmetric  steady 
solution. 

In  the  present  study,  we  use  a  free  slip  boundary  condition  on  the  western 
boundary.  As  Berio ff  and  Me  Williams  [3]  showed,  the  presence  of  a  no-slip  boundary 
condition  tends  to  destabilize  western  boundary  currents.  Eddies  produced  by  such 
currents  can  be  advected  towards  the  middle  of  the  domain  and  may  interact  with 
rings  that  propagate  towards  the  western  wall  from  the  model’s  interior.  Hence,  the 
nature  of  low  frequency  variability  and  of  multiple  dynamical  regimes  might  change 
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with  the  use  of  no-slip  boundary  condition  due  to  more  complicated  dynamics  of  the 
western  boundary  current.  This  problem,  however,  is  not  addressed  in  the  present 
thesis. 

In  Chapter  5,  we  describe  the  results  of  the  numerical  experiments  with  a  two- 
layer,  colliding  jets  model  which  shows  strong  dependence  of  the  flow  structure  on  the 
bottom  friction  parameter.  For  small  bottom  friction,  which  implies  a  physically  rel¬ 
evant  lower  layer  flow,  the  model  stays  in  a  unique  dynamical  regime  associated  with 
the  well-developed  recirculations  in  both  layers.  To  understand  the  nature  of  low  and 
high  frequency  variability  observed  in  the  model,  an  extended  spatial  EOF  analysis 
was  performed  for  filtered  streamfunction  fields.  The  low  frequency  spatial  EOFs 
are  found  to  explain  decadal  and  longer  variability  in  the  position  of  the  eastward 
jet  axis,  as  well  as  weak  wrapping  and  zonal  pulsations  of  the  recirculation  gyres. 
For  quasi-annual  and  mesoscale  variability,  the  EOF  analysis  shows  the  importance 
of  wave-like  structures  sitting  on  the  axis  of  the  jet  between  the  recirculating  gyres 
with  the  length  scale  of  1000  km  and  the  time  scale  of  8  to  9  months.  The  results  of 
the  EOF  analysis  for  the  regional  model  differ  from  the  stability  calculations  of  the 
zonally-symmetric  jet,  where  the  time  scales  are  much  shorter:  typically,  between  2 
and  4  months. 

For  large  values  of  the  bottom  friction  coefficient,  the  lower  layer  flow  becomes 
strongly  frictional,  and  the  recirculations  disappear.  However,  the  dynamics  of  the  up¬ 
per  layer  flow  becomes  more  complicated  displaying  both  multiple  dynamical  regimes 
and  low  frequency  variability  which  closely  resemble  those  of  the  reduced-gravity 
model  with  weak  interfacial  friction.  Such  similarity  between  the  reduced-gravity 
and  the  two-layer  models  is  possible  due  to  forcing  provided  by  the  lower  layer  flow 
on  the  upper  layer  flow,  associated  with  strong  bottom  friction. 

The  results  of  the  two-layer  calculations  show  that  for  physically  feasible  values 
of  the  bottom  friction  coefficients,  the  presence  of  baroclinic  instability  in  the  recir¬ 
culation  area  prohibits  the  existence  of  multiple  dynamical  regimes.  However,  weak 
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low  frequency  variability  associated  with  changes  in  the  structure  and  the  strength 
of  the  recirculating  gyres  is  still  possible.  In  contrast  with  the  reduced-gravity  case, 
where  steady  antisymmetric  solutions  are  useful  for  understanding  the  existence  and 
the  structure  of  a  high  energy  state  and  the  associated  dynamical  regime,  an  attempt 
to  find  steady  antisymmetric  solutions  by  employing  a  two-layer  model  with  enforced 
antisymmetry  failed.  In  disagreement  with  the  original  time  dependent  model,  where 
the  zonal  extent  of  the  recirculations  varies  between  500  km  and  2500  km,  the  recir¬ 
culations  reach  all  the  way  to  the  eastern  sponge  layer,  and  the  flow  is  characterized 
by  a  statistically  steady  nonzero  antisymmetric  eddy  field.  Therefore,  the  total  en¬ 
ergy  level  never  reaches  a  steady  state.  One  of  the  reasons  why  the  antisymmetric 
steady  state  does  not  exist  may  be  that  the  upper  layer  dissipation  mechanism  is 
presented  only  by  biharmonic  diffusion,  which  is  not  efficient  at  dissipating  energy. 
Therefore,  the  balance  between  the  nonlinear  terms  and  dissipation  is  not  possible, 
and  time  variability  is  essential.  The  question  still  remains  whether  the  correspond¬ 
ing  steady  states  are  relevant  to  the  behavior  of  the  time-dependent  model.  As  the 
system  becomes  progressively  more  complex,  as  in  the  case  with  both  barotropic  and 
baroclinic  instability  at  work,  the  model  might  become  purely  turbulent,  and  the 
resulting  dynamical  state  may  be  a  strange  attractor  for  a  given  parameter  regime. 

In  contrast,  in  the  two-layer,  wind-driven,  quasi-geostrophic  model  by  Berloff 
and  McWilliams  [4],  rare  transitions  between  three  distinct  dynamical  states  were 
possible.  Compared  with  a  corresponding  reduced-gravity,  wind-driven  model,  the 
presence  of  baroclinic  instability  weakened  and  regularized  low  frequency  variability. 

When  interpreting  the  results  of  the  reduced-gravity  and  the  two-layer  col¬ 
liding  jets  models,  one  must  remember  that  they  are  highly  idealized.  One  of  the 
basic  simplifications  used  in  the  present  study  is  the  quasi-geostrophic  nature  of  the 
flow,  which  overemphasizes  the  importance  of  horizontal  structure  of  the  ocean  and 
which  allows  for  only  small  perturbation  of  pressure  and  density  surfaces  from  their 
original  position.  And  although  the  quasi-geostrophic  approximation  proved  to  be  a 
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useful  tool  in  studying  the  mid-latitude  ocean  dynamics,  a  primitive  equation  model 
might  halt  low  frequency  variability.  Even  within  the  quasi-geostrophic  approxima¬ 
tion,  consideration  of  continuously  stratified  fluid  might  further  reduce  low  frequency 
variability. 

Low-frequency  variability  is  well  documented  in  observations  of  the  world 
oceans.  Twenty  two  years  of  hydrographic  data  along  137°£'  show  strong  inter-annual 
bimodal  path  variations  in  the  Kuroshio  (Qiu  and  Joyce  [30]).  The  net  transport  of 
the  Kuroshio  system  which  includes  the  jet  and  its  recirculations  increases  by  30% 
during  the  meander-path  years.  The  increase  is  primarily  due  to  the  decrease  in  the 
recirculation  strength.  Thus,  recirculations  are  strongest  during  the  straight-path 
years.  Observations  of  East  Australian  current  ( Roemmich  and  Cornuelle  [33])  and 
of  Brazil  and  Malvinas  currents  ( Olson  et  al.  [26])  also  display  strong  inter-annual 
variability.  Observations  of  sea  surface  height  in  the  Gulf  Stream  area  by  Kelly  et 
al.  [19]  show  a  trend  over  two  years  from  a  state  with  a  weakly-meandering  jet  and 
strong  recirculations  to  a  state  with  a  strongly  meandering  jet  and  small  recircula¬ 
tions.  Sea  surface  temperature  observations  by  Lee  and  Cornillon  [21]  indicate  annual 
and  inter-annual  oscillations  in  the  position  of  the  Gulf  Stream  path. 

Although  the  colliding  jets  model  isolating  the  internal  dynamics  of  the  jet  as 
sources  of  low  frequency  variability  is  highly  idealized,  it  produces  similar  variability 
for  a  wide  range  of  dissipation  parameters:  namely,  pulsations  in  the  zonal  extent 
of  the  recirculations,  weak  wrapping  of  the  recirculations  around  each  other  and 
meridional  shifts  in  the  position  of  the  jet.  Another  interesting  feature  of  the  model 
is  the  importance  of  wave-like  structures  with  time  scales  between  9  months  and  1.3 
years.  Lee  and  Cornillon  [21]  indicated  the  importance  of  a  9-month  cycle  in  the 
meandering  intensity  of  the  jet.  Clearly,  the  zonally  symmetric  model  was  not  able 
to  reproduce  this  feature,  since  most  energetic  wave  had  much  shorter  time  scales. 

Perhaps,  in  the  real  ocean,  there  are  several  mechanisms  that  could  be  re¬ 
sponsible  for  structural  changes  in  the  jet  and  in  the  recirculation  area,  such  as  the 
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presence  of  time  varying  winds,  complicated  topography,  the  instabilities  of  the  jet 
and  the  inherent  dynamics  of  the  steady  wind-driven  circulation.  The  colliding  jets 
model  isolates  the  importance  of  one  of  them,  the  internal  dynamics  of  the  jet,  and 
shows  the  differences  between  baroclinic  and  barotropic  instabilities  as  sources  of  low 
frequency  variability.  Further  studies  are  necessary  for  understanding  of  the  role  of 
time-varying  winds  and  topography  in  low-frequency  variability. 
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